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Abstract We address the problem of the exact computation of two joint spectral 
characteristics of a family of linear operators, the joint spectral radius (in short 
JSR) and the lower spectral radius (in short LSR), which are well-known different 
generalizations to a set of operators of the usual spectral radius of a linear operator. 
In this article we develop a method which - under suitable assumptions - allows to 
compute the JSR and the LSR of a finite family of matrices exactly. We remark 
that so far no algorithm was available in the literature to compute the LSR exactly. 

The paper presents necessary theoretical results on extremal norms (and on 
extremal antinorms) of linear operators, which constitute the basic tools of our 
procedures, and a detailed description of the corresponding algorithms for the 
computation of the JSR and LSR (the last one restricted to families sharing an 
invariant cone) . The algorithms are easily implemented and their descriptions are 
short. 

If the algorithms terminate in finite time, then they construct an extremal 
norm (in the JSR case) or antinorm (in the LSR case) and find their exact val- 
ues; otherwise they provide upper and lower bounds that both converge to the 
exact values. A theoretical criterion for termination in finite time is also derived. 
According to numerical experiments, the algorithm for the JSR finds the exact 
value for the vast majority of matrix families in dimensions < 20. For nonnegative 
matrices it works faster and finds JSR in dimensions of order 100 within a few 
iterations; the same is observed for the algorithm computing the LSR. To illustrate 
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the efficiency of the new method we are able to apply it in order to give answers 
to several conjectures which have been recently stated in combinatorics, number 
theory, and the theory of formal languages. 

Keywords Linear operator • joint spectral radius • lower spectral radius • 

algorithm • polytope • extremal norm • antinorm. 
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1 Introduction and background 

The joint spectral characteristics of linear operators are now applied in many 
areas, from functional analysis and dynamical systems to discrete mathematics 
and number theory. We focus on two characteristics: the joint spectral radius and 
the lower spectral radius, and elaborate a method of their exact computation 
applicable even for relatively high dimensions. 

The joint spectral radius of a set of matrices is a measure identifying the highest 
possible rate of growth of the norm of products of matrices (with no ordering 
and with repetition permitted) of the set. In contraposition, the lower spectral 
radius defines the lowest possible rate of growth. Both measures appear in several 
applications (see e.g. Strang [Str ). 

In this paper we consider the problem of the computation of both joint spectral 
characteristics for a finite set of matrices. In contrast to the fact that in the last 
twenty years much effort has been devoted to the computation of the joint spectral 
radius, very little is known about computing the lower spectral radius (to the best 
of our knowledge, the only available method of its approximate computation was 
presented in [PJB]). 

The joint spectral radius originated with Rota and Strang in 1960 [RSj. and 
became extremely popular after Daubechies and Lagarias ^DLi revealed its role 
in the study of refinement equations and wavelets. Since then it has found ap- 
plications in functional equations, approximation, probability, combinatorics, etc. 
(see jJl lPJB| for the extensive bibliography). Let 

M = {Ai,...,A^} 

be a finite family of linear operators acting in R*^. We write 

M'' = {Ad,... Ad, \d, e {l,...,m}, j = l,...,fc} 

for the set of all products of length k of operators from A4. The joint spectral 
radius (JSR) of the family A4 is 

p{M) = lim max ||B|| (1) 

fc— >oo BeM'' 

This limit exists for every family A4 and does not depend on the norm in R'* |BW| . 
Clearly, if A4 consists of one operator Ai, then p{A4) = p{Ai), where p{Ai) is 
the (usual) spectral radius of Ai, which is the maximal modulus of its eigenvalues. 
For any family A4 there is a positive constant ci such that 

max \\Ad, ■ ■ • AdJI > a p^ 
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for every fc £ N. The family is called non-defective if the inverse estimate holds, 
that is there is a constant C2 such that maxdj_..._d^ W^d^ ■ ■ ■ \\ < C2 p'' ■ It ap- 
pears that if a family A4 is irreducible, i.e., its operators do not share a common 
nontrivial invariant subspace of R'', then it is non-defective [PI] . Thus, for an 
irreducible family one has maxdi,...,tij, \\Ad^ ■ ■ .^t/JI ^ , where the symbol x 
denotes asymptotic equivalence ( x hk if there are constants ci , C2 > such 
that ciOfe < fefc < C2afc for all fc). Whence, the joint spectral radius is the expo- 
nent of polynomial growth for the largest norm of operator products of length k. 
The geometric sense of JSR is the following: p < 1 if and only if there exists a norm 
in such that || < 1 for all j = 1, . . . , m, where \\A\\ = sup||^||<;j^ is the 

corresponding operator norm. In other words, p < 1 precisely when there is a norm 
in M'^ with respect to which all operators from M are contractive. So, it is natural 
to expect that each family of operators possesses some special norms related to 
JSR. The following theorem established in 1988 by Barabanov [B] shows that this 
is indeed the case, at least for irreducible families. A norm || • || in M'^ is called in- 
t;aria7it for if there is a number A > such that maxj=i_..._m ll^ja^ll = -^||a^|| for 
every x G R''. It is shown easily that for every invariant norm one has A = p{M). 

Theorem 1 JS^ Every irreducible family M possesses an invariant norm. 

In practice it suffices to get a special norm with some weaker requirements, the 
so-called extremal norm. 

Definition 1 A norm \\ ■ \\ is called extremal for A4 if ||ylja;|| < p||a;|| for all 
X e R"*. 

Thus, the norm is extremal if and only if \\Aj\\ = p. Indeed, from 

the definition it follows that maxj=i_...^m \\Aj\\ < p; on the other hand, the sub- 
multiplicativity of operator norms yields maxj=i_...^rn || > P- Whence, this 
inequality becomes an equality precisely for extremal norms. This property justi- 
fies the term "extremal" . 

Clearly, any invariant norm is extremal, but not vice versa. Let us remark that 
for every extremal norm one has maxdj^...^(ij. || A^i^, . . . Aii-^ \\ = p^ , fc G N, i.e., the 
asymptotic equality becomes a sharp equality for all k. In particular, for fc = 1 
we have maxj \\Aj\\ = p. Thus, if we know an extremal norm, then we have the 
exact value of JSR. The main idea of the approach presented in this paper is to 
find JSR by constructing (in an iterative way) an extremal norm. 

We first present an algorithm for general sets of matrices, which under some 
suitable assumptions is able to check if a certain product in the multiplicative 
semigroup is spectrum maximizing. The algorithm is based on the computation of 
an extremal norm whose unit ball is a balanced polytope, and we provide it by a 
new criterion assuring a finite time termination and a new stopping condition. 

Then we analyze sets of matrices having an invariant cone (a most important 
case is given by families of nonnegative matrices). For such sets we refine the 
algorithm for the general case and exploit the invariant property of the set in 
order to make the algorithm faster. Under this assumption we are also able to 
determine an algorithm for the exact computation of the lower spectral radius, 
which appears to be the first algorithm able to provide an exact value of this 
important measure. 

The algorithms compute respectively a bounded and an unbounded polytope 
which represent the unit balls of, respectively, an extremal norm and antinorm for 
the considered set. 
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We write the formal routines of the algorithms and illustrate their efficiency by 
suitable examples and by numerical tests with randomly generated matrices. As we 
shall see, the algorithms find the exact value of JSR for general families of matrices 
(under some minor restrictions) in dimensions up to 20. For nonnegative matrices 
they work surprisingly fast even in dimension d = 100 and higher. Let us remark 
that our approach does not apply successfully to all families. There are cases, in 
fact, where the algorithms we propose are not able to finitely compute the exact 
values of the considered joint spectral characteristics, but only to approximate 
them. 

In view of negative complexity results for the problem of JSR computation [BT] , 
this is unlikely that there are effective methods applicable for all families of op- 
erators. Nevertheless, we claim that our approach works for the vast majority of 
families. The results of many numerical tests with randomly generated matrices 
of dimension from 5 to 100 (some of them are presented at the end of this paper) 
confirm this claim. In all the cases the algorithms found the exact values of JSR. 
As a further confirmation of this, we are able to apply the new method to solve 
several open problems in combinatorics and discrete mathematics. 

In the literature there are several methods for the computation of the JSR. 
Some of them work only for small dimensions d, but give either an exact or a very 
accurate value of p. For example, the method of polytope norms [Pl GZl BJPl 
IGZ2l[CGSZj : see also special methods in |VllP3llHMRj elaborated for particular 



Other methods aim to an approximate computation, such as the Kronecker 
lifting method [P2,BN , ellipsoidal norm method |BNTj . Gripenberg's branch- 
and-bound method [DLlfG] can work for bigger dimensions (mostly, up to 20), 
but produce pretty rough estimations. Recent approaches involving some modern 
tools of convex optimization (conic and semidefinite programming, sum-of-squares 
approximation, etc.) have rather good accuracy for higher dimensions (10 or even 
bigger) [PJIIPJB] . Most of those methods are actually based on the same simple 
idea. For each k we have 



The right hand side of this inequality converges to p as fc — > cxd, which follows 
from the definition. The upper limit of the left hand side also equals to p [BWj. 
So, choosing k large enough, it is possible to approximate JSR as close as we 
need. However, this possibility is purely theoretical, because in most of practical 
cases the number k grows as where e > is the relative accuracy of the JSR 
approximation, and C > is a constant, which may be large for high dimensions d. 
That is why the number of matrix products of A^*^ to look over becomes enormous. 
The reason is that the norm || • || in the right hand side of ([2]) may not suit our 
family A4, i.e., it may be far from the extremal norm of that family. That is why, 
to achieve a good approximation of JSR one needs to find an appropriate norm 
in R'^ for the right hand side of ([2]). Actually, all the methods of JSR computation 
use various techniques to find such a special norm for a given family A4. Those 
are, for instance, a polytope norm [Pl.GZlj, an ellipsoidal norm [BNTj . a norm 
generated by a cone [PJB] , a norm defined by a sum-of-squares polynomial [ PJ] , 
etc. Sometimes this idea leads even to finding the precise values of JSR. This 
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i/fc 
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happens when both the inequahties in ([JJ become equahties. If we write 

Pfe = max [p{B)Y^% 
BeM" , s<k '- ^ 

then for the extremal norm we have ||A|| < pk for aU A £ A4, and therefore 
Pk = P- For example, if all the matrices of A4 are symmetric, then the Euclidean 
norm is extremal, if they are all column-stochastic, then the Li-norm is so. In 
some practical cases people succeed in finding extremal norms for concrete pairs 
of matrices arising in various applications: Gripenberg [G^ (matrices of Daubechies 
wavelets, dimensions from 4 to 7), Hechler, Mofiner, and Reif |HMRj (matrices of 
the four-point subdivision schemes, d = 4), Protasov |P4j (de Rham matrices, 
d = 2), Guglielmi, Wirth, and Zennaro [GWZj (matrices of the Blondel-Theys- 
Vladimiov family, d = 2), Protasov |P3] (matrices of the binary partition function, 
d = 4, . . . , 12), Villemoes [V] (matrices of refinement equations, d = 2), Guglielmi, 
Manni and Vitale |GMVj (matrices in Hermite subdivision schemes) etc. The JSR 
computation in each case was a nontrivial problem and required special tricks 
applicable only for some narrow classes of matrices. 

The method of exact JSR computation presented in this paper is related to 
previous works (see [PlllGZlllBJPIIGZ2llCGSZj ') and aims to develop further ideas 
both for the general case, which we are goind to recall, and for certain specific 
important cases, like that of nonnegative matrices. The method is applicable for 
all families of matrices, under some general assumptions. The main idea proposed 
in the above mentioned papers is to build an extremal norm, whose unit sphere 
is a polytope. At the first step we look over all products of matrices from M of 
length at most I, and find a product /7, for which the value [p{n)Y^^ is maximal 
(n is the length of 77). Then we denote this value by pi and try to prove that 
p{M) = PI. 

Definition 2 A product U 6 A4^ is a spectrum maximizing product (s.m.p.) if 

[p(77)]^/" = piM). 

To prove that 77 is an s.m.p. it suffices to have an extremal norm || • || in R'^, 
for which \\Aj\\ < pi , j = I, . . . ,m. By Q in this case we indeed have pi = p. 
We try to build a polytope extremal norm, whose unit sphere is some polytope P. 
Such a polytope will also be called extremal. It is characterized by the property 
AjP C piP , j = l,...,m. The polytope is constructed successively: its first 
vertices are the leading eigenvector vi of 77 (i.e., the eigenvector corresponding to 
the largest by modulo eigenvalue, which is assumed to be real for the moment), 
the leading eigenvectors Vi of the {n — 1) cyclic permutations of 77, and the same 
vectors taken with minus, i.e., —Vi. We call an eigenvalue A of an operator A 
leading if |A| = p{A). 

Then we consider their images {pi)~^ AjVi , j = 1, . . . ,m and remove those are 
in the convex hull of the previous ones, etc., until we obtain a set of points V such 
that 

{pi)-^A,V C co4V), j = l,...,m. 

By cOs(V) we denote the symmetrized convex hull: cOs (V) = co (VU (—V)), where 
co(-) is the (usual) convex hull. Then the polytope P = cOs(V) possesses the 
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desired property: {pi)~^AjP C P, so P is an extremal polytope. This implies 
p = pi. The algorithm involves standard tools of linear programming. 

In case the leading eigenvalue of 77 is complex, one has to replace polytopes 
by the so-called complex polytopes (see e.g. |GZ3] ). 

Our goal is to develop this approach for general families of matrices in higher 
dimensions, to analyze the structure of extremal polytopes and to derive the con- 
ditions of convergence of this algorithm. 

Let us now emphasize the shortcomings of our approach. First of all, not ev- 
ery family of matrices has an s.m.p. Moreover, even if a non-defective family A4 
possesses an s.m.p., it may not have extremal polytopes (neither real nor com- 
plex ^JPj). For such families our method apparently does not work. Another dis- 
advantage appears, when the s.m.p. is not unique, up to cyclic permutations. In 
this case an extremal polytope, even if it exists, in general cannot be found by 
our method. The first two cases are rather pathological. It required constructing 
special nontrivial examples to show that they are possible |BTVI[JP] . The third 
case of multiple s.m.p., in contrast, being also quite rare in general, nevertheless, 
appears in practical applications. 

We believe that our method can be extended to this case as well, which may 
be a challenging problem for further research. 

To work with those "bad cases" , we apply our approach also to approximate 
computation of JSR. The algorithm constructs a polytope, which is either extremal 
or not. If it is, then the JSR is found. Otherwise, we stop the algorithm after a 
certain iteration, say the A''-th, and use the obtained polytope as a unit ball of 
the corresponding norm in estimations ([2]). In most cases this gives very sharp 
bounds for JSR. Thus, for an arbitrary family A4 the algorithm either produces 
an extremal polytope, or a polytope norm that gives good upper and lower bounds 
for JSR. Proposition [T] guarantees that both those bounds converge to p{A4) as 
N ^ oo. 

The second part of the paper deals with the lower spectral radius (LSR) defined 
as follows: 

p{M) = lim min ||73||i/\ (3) 

fc— >cxj BEM'' 

Thus, LSR is the exponent of asymptotic growth of the minimal product of opera- 
tors from the family A4. This notion defined in [Gu] have been applied in problems 
of dynamical systems, functional analysis, coding theory, combinatorics, number 
theory, etc. (see ^ for many references). The limit in ^ always exists and does 
not depend on the norm. A simple observation is that LSR can be estimated by 
the usual spectral radii as follows: 

p{M) < min [p(S)l'/'' < min ||B|| (4) 

In contrast to inequality Q for JSR, estimation gives only upper bounds. 
In fact, there is no effective lower bounds for LSR, and this causes the main 
difficulty for its computation. Basically, the lower spectral radius is still harder 
to compute or to estimate than the joint spectral radius (see, for instance, [TBj 
for the corresponding complexity results). The notions of invariant and extremal 
norms cannot be directly extended to LSR. The reason is that the operation of 
taking minimum of several functions, in contrast to the maximum, does not obey 
convexity. This means that the pointwise minimum of several convex functions 
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may not be convex. Hence, the function f{x) = minj=i,...,m , a; G M , in 

general, is not a norm in M''. To overcome this difficulty, we use in Section|6]a notion 
of antinorm defined on a convex cone K (ZW^ (Definition|31). This notion originated 
in jP6j to study the Lyapunov exponents of linear operators. As we shall see, it can 
also be applied to analyze the lower spectral radius. We prove that every family of 
operators that share a common invariant cone K possesses an extremal antinorm 
on that cone (Theorem [5]). This allows us to extend the new approach to the LSR 
computation, replacing norms by antinorms, and polytopes by infinite polytopes, 
i.e. the sets of the type co (V) + K, where V C R'' is a finite set, and C M'^ is 
a cone. In particular, this approach can be used for nonnegative matrices, since 
the corresponding operators preserve the cone K = W^. . This yields an algorithm 
of exact computation of LSR for nonnegative matrices. In numerical examples we 
show that the algorithm works well for rather big dimensions (like d = 100). Let us 
note that the problem of LSR computation for nonnegative matrices arise naturally 
in combinatorics, discrete mathematics, and number theory [dlP3l[TPB2] . Some 
of those applications will be considered in detail in Section [9l See also [MSllFVj 
for applications to the problem of stabilization of switched linear systems. 

The main results of the paper can be summarized as follows: 

(i) we analyze the considered algorithm for the JSR computation of an arbitrary 
family and improve it by elaborating a stopping criterion that indicates whether 
a chosen product 77 can be an s.m.p. or not. If our initial guess is wrong, and 77 
is not an s.m.p., then the criterion determines it (usually, after a few iterations) 
and suggests a new candidate for s.m.p. with a bigger spectral radius. 

(ii) Theorem |3] in Section [5] gives a criterion for a family M insuring that the 
algorithm terminates within finite time, i.e., produces an extremal polytope. 

(iii) we improve the considered algorithm when applied to nonnegative matrices; 
the new algorithm finds the exact values of JSR in much higher dimensions 
(up to d = 100); 

(iv) we obtain a new algorithm which is able to exactly compute the LSR for 
families of nonnegative matrices, by computing a polytope extremal antinorm; 

(v) as examples we compute the exact values of JSR for special families of matri- 
ces (of dimensions up to 40) from well-known problems of combinatorics and 
number theory. This, in particular, allows us to solve three open problems. We 
discuss this aspect below in more detail. 

(vi) we provide numerical tests with randomly generated matrices (both arbitrary 
and nonnegative) , showing that for all considered cases the algorithms produce 
extremal polytopes and, consequently, the exact value of the JSR (LSR). 

The structure of the paper is the following. We describe the algorithm for 
JSR computation in three possible cases, which will be considered separately and 
called (R), (C) and (P). The case (R), when the leading eigenvalue of the prod- 
uct 77 e (a candidate for s.m.p.) is real, is recalled and further analyzed in 
Section [2] We discuss an algorithm for constructing an extremal polytope and for 
computing JSR, give necessary explanations and proofs, and establish two effi- 
ciency results: on the stopping criterion (to indicate within finite time, whether 
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the chosen product 77 is s.m.p. or not) and on the estimation for JSR. Thus, Al- 
gorithm (R) either terminates within finite time, in which case JSR is found, or 
produces lower and upper bounds converging to JSR. According to our numerical 
experiments with randomly generated matrices (Section [9|), for almost all matrix 
families Algorithm (R) finds the exact value of JSR, and works efficiently for 
dimensions up to 20. 

In Section [3] we briefly consider the case (C), when the leading eigenvalue 
of 77 is complex, where we refer to |GWZ|[GZ2IIGZ3] . The algorithm and all the 
efficiency results are very similar, but with complex polytopes. By the numerical 
results in Section |9l it works slower than Algorithm (R) , and works in smaller 
dimensions. 

In Section[3]we consider the case (P), when all matrices are nonnegative. In this 
case the corresponding Algorithm (P) works faster and much more efficiently. Of 
course, (P) is a special case of (R), which is, in turn, a special case of (C). In fact 
all the three algorithms are very similar and differ in a few key details. Nevertheless, 
we describe them separately and independently of each other for convenience of 
the reader. Besides, their practical efficiency is very different, and it would be 
non-reasonable to compute JSR of nonnegative matrices by Algorithm (R) or 
by (C). 

In Section[5]we formulate one of the main results of the paper. This is a criterion 
of terminating of Algorithms (R), (C), and (P) within finite time (Theorem |4}. 
It shows that an algorithm produces and extremal polytope and finds the precise 
values of JSR if and only if the family A4 has a dominant product (Definition [3]) . 
In particular, if the algorithm terminates within finite time, then 77 is a dominant 
product for A4. 

In Section[6]we extend our method to the lower spectral radius computation. To 
this end we first define an antinorm, prove several theoretical results about it, and 
then describe Algorithm (L) for the exact computation of LSR of nonnegative 
matrices. Its practical efficiency for randomly generated matrices (Section [9|) is 
approximately the same as for Algorithm (P). 

Section [3 presents two detailed examples in dimension 2 to illustrate the algo- 
rithms. 

In Section [8] we consider applications to several problems of combinatorics, 
coding theory and number theory. In the problem of asymptotic growth of the 
number of overlap- free words f i)8.ip we compute precise values of exponents of the 
upper and lower growth. This proves two conjectures stated in 2008 |JPB2j . Then 
in !j8.2l we do the same for the problem of density of ones in the Pascal rhombus, 
and disprove one previously know conjecture. In and ?)8.4I we find precise 

values of the lower and upper growths of the Euler partition functions for some 
values of the parameters. 

Section[9] presents the results of numerical tests for JSR and LSR computation 
for randomly generated matrices of dimensions from 5 to 100. In all the cases 
the algorithms find the exact values of JSR and LSR, which suggests that our 
approach generically has finite convergence. 

In the sequel we assume that the basis {ei}f^i of the space R"^ is fixed and do 
not distinguish between operators and the corresponding matrices. An eigenvalue A 
is simple, if it is of multiplicity 1. The largest by modulo eigenvalue of an operator 
B is called leading and denoted by Amax (if there are several such eigenvalues, then 
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each of them is leading). We use the foUowing notation: B* is the operator adjoint 
to B, int M is the interior of a set M, co (M) is the convex hull of M. We use the 
short abbreviation "LP" for linear programming problems. 

2 Computing of the joint spectral radius: the case of real leading 
eigenvectors (R) 

In this section we present Algorithm (R) for JSR computation. 

We consider an irreducible family A4 = {Ai, . . . , Am}- For some (as large 
as possible) I we look over all products 77 of length < I and take one with the 
biggest value [p(77)]^/", where n is the length of the product. We denote it as 77 = 

Ad,, • • • 4^1 • ^ ^ 

Let M = {Ai, . . . , Am} be the normalized family, where Ai = [/o(77)] At. 
For the product 77 = Ad^ ■ ■ ■ Ad-^ we have p(77) = 1 which implies p{A4) > 1. 

Define, for an arbitrary nonzero vector d € M'' the set 

n{v)= U {rv \reM^}, (5) 

fc>0 

(where = Id, the identity matrix), i.e. the set obtained by joining v to all vec- 
tors obtained by applying the products of the semigroup of A4 to v. The following 
theorem (see |P1| and |GZip relates the set f2(w) and an extremal norm for Ai. 

Theorem 2 Let M = {Ai, . . . , Am} be irreducible and such that p{A4) > 1 
and let i7(w) (for a given v ^ 0) be a bounded subset of M."^ spanning M."^ . Then 
p{A4) = 1. Furthermore the set 

CO, {n{v)) = co{n{v) U -n{v)) (6) 
is the unit ball of an extremal norm \\ ■ \\ for AI (and for A4 ). 

The main idea of the algorithm we present is to finitely compute the set ([6]) 
whenever it is a polytope. Let us clarify this key point. 

We say that a bounded set 7-" C R'' is a balanced real polytope (b.r.p.) if there 
exists a finite set of vectors V = {ui}i<i<p (with p > d) such that span(V) = M*^ 
and 

P = co,{V) = co{V,-V). (7) 

Therefore 

7-" = = ^ tx X with — qx < tx < Qx, Qx > O^x eV and ^ (Jx < l|. 

xev xev 

The set P is the unit ball of a norm || • ||p on R'^, which we call a real polytope 
norm. 

Assume that the hypotheses of Theorem [2] hold. The possibility of actually 
determining an extremal polytope norm, if any, crucially relies on the search of the 
initial vector v, which we will address later in Theorem 2] that suggests to choose 
u as a leading eigenvector of 77 (although a different choice would be admissible). 
We will assume here v to be real. 
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The idea is that of computing the set Q{v) by applying recursively the family 
M to a finite set of vectors (which in the beginning is simply the vector v), checking 
at every iteration h whether M maps the symmetrized convex hull {cOs{n^-'^{v))) 
of the computed set of vectors 

n''-\v)= U {rv \reM''j, 

0<k<h-l 

into itself. 

Algorithm (R) we are going to present is similar to the one described in [GZl] ; 
the main differences are that a new vertex is included even if it lies on the boundary 
of the current polytope, all the leading eigenvectors are considered as starting ver- 
tices of the searched extremal polytope and a new and efficient stopping criterion 
is added. 

We start by an auxiliary result and then describe the algorithm. 

Lemma 1 Let an operator B have a unique simple leading eigenvalue A G M with 
the leading eigenvector v; let also v* be the leading eigenvector of B* such that 
{v*,v) = 1. If for some operator C one has |(?j*,Cw)| > 1, then for sufficiently 
large r the operator B^C has a unique simple leading eigenvalue, which is real and 
bigger than A by modulo. 

Proof. Without loss of generality, after a suitable normalization, it can be assumed 
that A = 1. Since all other eigenvalues of B are smaller by modulo than 1, it follows 
that B^ converges to the one-rank operator Soo x = {v*,x)v as r — ^ oo. Hence 
B^C converges to the operator BooC, whose unique simple leading eigenvalue is 
{v* , Cv), which exceeds 1 by modulo. □ 

Remark 1 If v and v* are the leading eigenvectors of B and B* respectively, then 
these vectors cannot be orthogonal, otherwise the leading eigenvalue is not simple. 
So, {v* ,v) ^ 0, and hence, after a suitable normalization it can always be assumed 
that {v* ,v) = 1. 

It is well-known that the problem of JSR computation has to be considered 
only for irreducible families of matrices, which do not possess common invariant 
linear subspaces. Otherwise this family is factorable in a suitable basis in R'^: 
all the matrices Aj get a block upper-triangular form, and p{A4) equals to the 
maximal JSR of the blocks. This reduces the problem of JSR computation to 
several problems in smaller dimensions. Therefore, in the sequel of this section we 
assume that M is irreducible. 



2.1 Algorithm (R) 

Initialization. Given the irreducible family M = {Ai, . . . , Am} we look over all 
products 77 of length < I and consider the shortest product 77 such that [/3(77)]^''" 
is maximal, where n is the length of the product. We denote it as 77 = Ad„ ■ ■ ■ Ad-^ 
and consider the main assumption: 

(i) The product 77 has a real nonzero leading eigenvalue. 
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We assume that the leading eigenvalue A max of n^is positive; the case of negative 
eigenvalue is considered in the same way. Let M = {Ai, . . . , Am} be the nor- 
malized family, where Ai = [p{n)] "^Z" Ai. For the product 77 = A^^ ■ ■ ■ A^-^ we 
have Amax = 1- 

Let 77i = 77 , 77, = • ■ ■ A^^A^^ ■ ■ ■ /l^. be a cyclic permutation of 77i, 

i = 2, . . . , n. We denote by vi the leading eigenvector of 77i, for which 77iwi = v\. 
If it is not unique (in which case A max is multiple) we take any of them. Then for 
every i > 2 we set 

Vi = Ad.._-^ ■ ■ ■ Ad^vi . 

Thus, Vi is a leading eigenvector of 77i. 

In case 77 has a unique simple eigenvalue, we also need the corresponding 
dual system of vectors: vl the leading eigenvector of the conjugate operator 77* 
normalized by the condition {vl ,vi) = l (see Remark [1} , and 

for i = 2, ... ,71. Thus, v* is the leading eigenvector of 77^, and {v*,Vi) = 1. If 
the leading eigenvalue of 77 is multiple or not unique, then we do not need the 
conjugate system. 

Set k = 0. We set Vo = Uq = {vi, . . . ,Vn} and TZq = {{vi, Ap)\ i = 
l,...,n; p= l,...,m, p^dr}. 

Main loop 

For k > 1. We have finite sets Vk-i C R''',llk-i C Vk-i, and TZk-i C 
Uk-ixM-SetVk = Vk-i And Uk = 0. Take an arbitrary pair («, A) € 7?.fc_i and 
compute the norm whose unit ball is the polytope with vertices Vk-i and —Vk-i, of 
the corresponding vector z = Av. This is done by solving the following LP problem 
with variables {tx}xeVk i {Qx}xeVk s-nd to (which represents the reciprocal of the 
value of the norm) : 

max to 

subject to to 2 = ^x 

xeVk 

<tx< > Vx e Vfc 

and J2 <l, 

The value of the problem, i.e., the value max to will be denoted by t^^ ^y. Thus, 
for a given pair {v,A) 6 TZk-i we obtain the value t^-^ ^y. 

If t^^ ^y > I, then we leave the sets Vk and Uk as they are, take the next 
pair {v,A) £ TZk-i and consider problem ^ for it. 

If t^^ < 1, then we distinguish between two cases 

If the leading eigenvalue of 77 is unique and simple, we apply the following 
Stopping criterion: 
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For a given pair {v, A) we check the condition 

\{v*,Av)\ < 1, j = l,...,n. (9) 

If ^ is satisfied, then we set Vk = { Av} ,Uk = UkU {Av}, take the 
next pair {v,A) £ TZk-i and consider problem ([8]) for it. 

Otherwise If ([9]) is not satisfied, then 77 is not an s.m.p. for A4, and 
p{M) > 1 (Lemnia[l]). We stop the algorithm and go either to the Final step, or 
back to the Initialization. In the latter case we need to find another candidate 
s.m.p. The first option is to increase / and to look over all products of a bigger 
length. Lemma [T] provides also a different approach. We take an index j, for which 
I {v* , Av)\ > 1. Applying Lemma [1] for the vectors v* and — v* , we conclude 
that there is r such that Amax(77j' As^ ■ ■ ■ Ag-^) > 1, where As^ ■ ■ ■ As-^vj = Av. 
We take the new initial product 77 = 77^ As^ ■ ■ ■ As-^ and restart the algorithm. 

End If 

Otherwise If the leading eigenvalue of 77 is not unique or multiple, then we 
do not apply the stopping criterion, and set Vk = VkU { Av} ,b(k = UkU {Av}, 
take the next pair {v,A) £ TZk-i and consider problem (|8]) for it. 

End If 

The kth step is over when the whole set TZk-i is exhausted. 

If Uk = 0, then p{M) = 1, and so p{M) = [p(77)]^/". The extremal polytope 
is Pk-i = cOs (Vfe_i), and the s.m.p. for A4 is 77. The algorithm terminates having 
performed k steps. 

Otherwise If Uk 0, then we set TZk = Uk x and continue. 

End If 

End For 

Final step. If the algorithm has not terminated, then we stop it after some 
N steps, denote tjv = jnin X\' where t, is the solution of LP 

problem ([5} for the last step, i.e., for k = N, and have the following estimate for 
the joint spectral radius of the family A4: 

[pW]'/" < p{M) < t],'[pm'/-. (10) 



End of Algorithm (R). 

Remark 2 An important difference with respect to previous similar algorithms is 
that if at step k a new vector v lies on the boundary of the polytope Pk-i = 
COs (Vk-i) (which means t^^ a } ~ ^ ^'^^ some p) then we include the vector as 
a new vertex. Clearly this condition is non generic and requires - to be tested in 
floating point arithmetics - the use of a suitable error tolerance. 

Before we give the proofs, let us explain the general scheme of the algorithm. 
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2.2 The cyclic tree structure of the algorithm 

Consider a combinatorial cyclic tree T defined as follows. The root is formed by 
a cycle B of n nodes vi, . . . ,Vn- They are, by definition, the nodes of zero level. 
For every i < n an edge (all edges are directed) goes from Vi to t^i+i, where we set 
Vi+i = vi. At each node of the root m — 1 edges start to nodes of the first level. 
So, there are n{m — 1) different nodes on the first level. The sequel is by induction: 
there are n{m — l)m''~^ nodes of the fcth level, fc > 1, from each of them m edges 
("children") go to m different nodes of the (fc + l)st level. 

Consider now an arbitrary word b = d„ . . . di of length n > 1, where each dj 
belongs to the alphabet {1, . . . , m}. The product of several words is their concate- 
nation. We assume that b is irreducible, i.e., is not a power of a shorter word. 
To every edge of the tree T we associate a letter d as follows: the edge viVi+i 
corresponds to , i = 1, . . . , n; at each node m edges start associated to m dif- 
ferent letters. To a given word ■ ■ ■ qi we associate the node, which is the end 
of the path from vi along the edges qi,q2, . . . ,qk- For example, the empty word 
corresponds to vi, the word b also corresponds to vi, the word d2di corresponds 
to V3, the word d2 corresponds to either V2, if d2 = di, or to a child of vi from the 
first level, otherwise. This tree is said to be generated by the word b, or by the 
cycle B. 

For a family of operators A4 = {Ai, . . . , Am} and for some product 77 = 
Ad„ ■ ■ ■ Ad-i with an eigenvalue 1 we associate the cyclic tree T generated by the 
word dn . . . di. The node vi corresponds to an eigenvector with the eigenvalue 1; 
to a given node v (zT we associate a point Aq^. . . . Aq-^vi, where the word qk ■ ■ ■ qi 
corresponds to the node vi. 

When we start the algorithm, we take the set B = {vi, . . . ,Vn} as the root of 
the tree. At the first step we take any node Vi and consider successively its (m — 1) 
children from the first level. For each neighbor u = Avi, where A £ M \ {^dj} we 
solve LP problem ^ and determine, whether u belongs to the interior of the set 
cOs (Vi), where cOs(M) = co{M, — M} is the symmetrized convex hull. If it does, 
then u is a "dead leaf" generating a "dead branch" : we will never come back to u, 
nor to nodes of the branch starting at u (so, this branch is cut off). If it does not, 
then u is an "alive leaf" , and we add this element u to the set Vi and to the set IA\ . 
After the first step all alive leaves of the first level form the set Ui . At the second 
step we deal with the leaves from lAi only and obtain the next set of alive leaves 
of the second level IA2, etc. Thus, after the fcth step we have a family Z//^ of alive 
leaves from the fcth level, and a set Vfc, = Uj^qUj. A node u belongs to Vk iff its 
level does not exceed fc and it belongs to an alive branch starting from the root. 
The polytope is the symmetrized convex hull cOs (Vfc). The polytope Pk-i is 
extremal iff Wfc = 0, i.e., the fcth step produces no alive leaves (only dead ones). 
This means that there are no alive paths of length fc from the root. Therefore 
Pk = Pk-i- Otherwise, \ilAk is nonempty, we make the next step and go to the 
(fc -h l)st level: take children of each element of Uk, determine whether they are 
alive or dead and proceed. 
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2.3 Explanations and proofs 

The algorithm produces a sequence of embedded polytopes fi C -P2 C . . . such 
that Pj+i = cOs {AiPj , . . . , AmPj } for every j. If the algorithm terminates after 
the fcth step, then Pk = Pk-i- The fcth step is actually needed only to ensure 
that the polytope Pk-i is extremal, i.e., Aj Pk-i C Pk-i , j = 1, . . . ,m. In this 
case Pk-i possesses an interior of nonzero measure, otherwise its linear span is a 
common invariant nontrivial subspace of the family M, which contradicts the irre- 
ducibility assumption. Moreover, Pk-i is centrally-symmetric, hence G intPfc_i. 
This, in particular, yields that if for some £ R'^ and t > 1 one has tv £ Pk-i, 
then V £ int_Pfc_i. Thus, if the value t of LP problem ([8]) is bigger than 1, 

then Av 6 intPfc_i. Thus all dead leaves removed by the algorithm are internal 
points for Pk-i- 

In the Minkowski norm || • ||fc_i whose unit ball is given by Pk-i, one has 
||A||fc_i < 1 for aU A e M, therefore p{M) < 1. On the other hand, p(M) > 
p{ny/" = 1, hence p{M) = p{ny/'' = 1, an so p{M) = piny/"". Thus, if 
the algorithm terminates within finite time, then the s.m.p. and the exact value 
of JSR are found. 

Suppose the algorithm does not terminate within finitely many steps. After 

the final step we take the polytope Pn-i = cOs (Vjv-i) as a unit ball of the new 

norm 11 ■ ||jv-i in R"*- Then max ||A||jv-i > p(M). We have 

AeM 

max||^|U_i = [p(77)i/"] max = [p(77)i/"] • max (t^^ 

where t^^ is the value of LP problem ^ for k = N . Therefore, p[M.) < (tiv)~^, 
and after multiplying by [p(77)]^^" we arrive at pO|) . 

Remark 3 Although we perform operations numerically, the obtained results have 
to be considered exact since apart from the vertices of Pk-i, all other vectors 
obtained by applying the scaled matrices Aj to the vertices are either vertices or 
internal points to the polytope Pk-i- 

Remark ^ By the construction of the algorithm, each vertex of the polytope Pk 
belongs either to Vk or to — Vfc. However, not all elements of the set Vk are actually 
vertices: some of them may lie in the convex hull of the others. This means that 
in general the set Vk is not an essential system of vertices. Nevertheless, for the 
sake of simplicity we call all elements of Vk vertices. 

Remark 5 Actually Algorithm (R) can be applied to a reducible family Ad as well. 
If the algorithm terminates after fcth iteration, and the set Vk C US'* does not lie 
in a linear subspace of a smaller dimension (i.e., the system of equations {x,v) = 
, ?J G Vfc has only trivial solution a; = 0), then Pk-i is an extremal polytope, and 
77 is an s.m.p. Thus, one can apply Algorithm (R) without preliminary checking 
of irreducibility of M. Nevertheless, if the family Ad is reducible, then it is always 
better to factorize Ad before applying Algorithm (R), because this reduces the 
dimension of matrices. 
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2.4 Efficiency resuits for Aigoritlim (R) 

If tlie aigoritlim terminates, then it finds the exact value of JSR, otherwise esti- 
mate (|IOp gives its approximate value with the relative error e = (ijv)~^ ~ 1- 
This error depends on two integer parameters: the maximal length I of the prod- 
ucts, among which we choose an s.m.p. 77, and the number of iterations N of the 
algorithm. Let us show that the error s tends to zero as both these parameters 
increase: 

Proposition 1 For an arbitrary irreducible family A4 we have tjv — > 1 as Z —> oo 
and N ^ oo. 

Thus, both sides of inequality (|I0|) tend to p{A4) as I — )■ oo , Af — > oo. Algo- 
rithm (R) either finds the value of JSR or provides lower and upper bounds for it; 
those bounds are arbitrarily close to each other, whenever both / and A'' are large 
enough. 

In the proof we use Dini's theorem on monotone convergence: if a sequence of 
continuous real- valued functions defined on a compact metric space Q is monotone 
and converges pointwise to a continuous function, then this convergence is uniform 
on Q (see Ru, theorem 7.13]). We use the Minkowski norm || ■ || d associated to a 
given symmetric convex body 7) C R*^ as follows: || • || d = inf {t~^ \ t > 0, tx £ 

Proof of Proposition [l] Assume first that p{A4) = 1, i.e., that 77 is an s.m.p. 
The algorithm produces the polytopes {PfcjfceN such that 

Pk C Pk+i = co^ {liPfc,...,1^7'fc}. 

Since the family A4 is irreducible, there is p > 1 such that all the polytopes Pk 
have nonempty interior for k > p. Hence, for k > p the polytope generates 
the Minkowski norm = II • II p,. For each x eR'' the sequence {fk{x)}k>p is 

non-increasing. Moreover, it is uniformly bounded below by a positive constarit^ 
because all the polytopes {Pk}k>p are contained in some ball, since the family A4 
is non-defective. Therefore, the sequence fk{x) converges pointwise to a function 
f{x), which is also a norm in R*^. By Dini's theorem, this convergence is uniform 
on any compact subset of R"^. In, particular, it is on the unit sphere S = {x G 
R'^ I f{x) = 1} of the norm /. Thus, fk{x) — > 1 uniformly for a; G S', as fc — > cxd. 
Hence, there is such that fN-i{x) < 1 + £ for all x £ S, whenever N > 
Ne. Consequently, (tw)"^ = sup^gp^ /Ar_i(a;) < sup^g ^ /at-i (x) < 1 -|- e, 
which completes the proof for the case p = I. Consider now the general case. 
We have [p(77)]^^" — > p{A4) as I — > cxd, where, let us remember, n = n{l) is the 
length of 77. Hence, for every S > there is Is such that p{A4) < I + 5. Since 
each polytope P^ continuously depends on the family A4 = [p{n)]~^^'^A4, for all 
sufficiently small 6 one has (tjv^)"^ < 1 + e. This inequality holds for all I > Is- 
It remains to note that for every family A4 the value tN is non-decreasing in N. 
Indeed, tj^Pj^ C 7-V-i, hence tj^AjPj^ C AjP^-i for every j = 1, . . . , m, and so 
tNPN+i C Pn- Therefore, tpf^i = sup{t > | £7-V+i C Pn} > ijv- We see that 
{tN)~^ < I + s, provided A'' > A^e- Thus, for every e > there are Is and A^e such 
that (^at)"^ — 1 < £, whenever N > and I > Is- 

□ 

Let us now show the efficiency of the stopping criterion. 
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Proposition 2 Assume the leading eigenvalue of U is real, unique and simple. If 
the assumption of the algorithm is wrong (i.e., II is not an s.m.p.) then for every 
j = 1, . . . ,n condition ^ is violated at some step. Conversely, if condition {PP is 
violated at some step for some j , then U is not an s.m.p. 

Proof . The sufficiency follows from Lemma [TJ To prove the necessity suppose 
p{M) > 1; then for every point w ^ and for every number R > there is a 
product C of operators of the family M such that ||Ct;|| > R [Pll theorem 1]. On 
the other hand, since the family M is irreducible, it follows thaVthere is 7 > such 
that for every point y ^ the set Kd{y) = cOs {Ay | A G A?'^} contains a ball 
of radius ^\\y\\ (see [K2| for the proof), where 7 > is a constant. Applying these 
results to the points v = vi and y = Cvi and using the fact that the polytope Ps+d 
contains Kd{y), where s is the length of the product C, we see that the polytope 
Ps+d contains a ball of radius ^R centered at the origin. Therefore, 

sup \{v*,x)\ > 7ii||?;*||. 

Since this supremum is attained at some vertex of Ps+d, which is produced by 
the algorithm, we see that condition ([9]) will fail by the (s + d)-th step, whenever 
i?> 1/(711^^,11*). 

□ 

Thus, if the chosen product 77 is not an s.m.p., then the stopping criterion 
always determines this in finite time. In Section [5] we formulate Theorem [3] that 
gives a sharp criterion for the algorithm to terminate in finitely many steps (and, 
respectively, to produce an extremal polytope) . 

3 Computing of the joint spectral radius: the case of complex leading 
eigenvectors (C) 

For the theoretical results and the algorithms relevant to this case we mainly 
address the reader to the papers |GWZl[GZ2l[GZ3] . 

We recall from GZ3, VZ^ the definition of a balanced complex polytope, which 
generalizes to the complex case a centrally symmetric real polytope. 

Let V = {i'i}i<i<p be a finite set of vectors, then 

absco(V) = |z e I 2 = ^ with ^ < l|. (11) 

xev xev 

Definition 1 A set P C is a balanced complex polytope (b.c.p.) if there exists 
a finite set of vectors V = {i'i}i<i<p such that 

span(V)=C'' and P = absco(V). (12) 

Moreover, if absco(V') £ absco(V) for all V' 5 V, we say that V is an essential 
system of vertices for V. Every vector uvi with u £ C, \u\ = 1, is called a vertex 
of -p. 

Note that geometrically a b.c.p. V is not a classical polytope (see |GZ3] V 
A polytope norm can be defined in a natural way. 
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Lemma 2 Any b.c.p. V is the unit ball of a norm \\ ■ Wv on . 
The proof is immediate (see e.g. [GZ3j l 

Definition 2 We shall call complex polytope norm any norm || • Hp whose unit 
ball is a b.c.p. V. 

The corresponding vector norm is characterized by the following Lemma (for a 
proof see [GZ3 ). 

Lemma 3 Let V be a b.c.p. and let \\ ■ W-p be the corresponding complex polytope 
norm. Then, for any z G C^, it holds that 

||2:||-p = I maxtp I ipz = t^ X, ^ Iti^l < l|, (13) 
xev x<ev 

where V = {wi}i<i<p is an essential system of vertices for V . 

Complex polytope norms are dense in the set of all norms defined on and 
consequently the corresponding set of induced matrix complex polytope norms is 
dense in the set of all induced d x d- matrix norms (see ^CZSj). This implies the 
following important property: 

p{M) = inf max \\A\\r 

where || • Hp denotes the set of polytope norms. 

From an algorithmic point of view, the above property has the consequence 
that although an extremal polytope norm may not exist, it is possible to compute 
a polytope norm which is e-close to an extremal one, for any e > 0. 



3.1 Main differences between Algorithm (C) and Algorithm (R). 

Algorithm (R) extends to the complex case in a direct way, as well as the con- 
vergence and approximation results (see |GWZ] . [GZ2]). The only (important) 
difference lies in the computation of the polytope norm of a vector. We obtain this 
by rewriting psp as a real optimization problem. 

Let V = absco(V) (with V = {vi,V2, . . . ,Vp}) be a b.c.p. and || • \\-p the asso- 
ciated norm. For any z £ C^, we write (|13p (with t^ = ax + il3x) in the following 
way: 

max to 

subject to Yl '^x Re(a:) — Im(a;) = toRe(z) 

J2 axIm{x) + PxRe{x) = toIm{z) ^^^^ 
xev 

and J2 PI < 1 

This problem can be efficiently solved in the framework of the conic quadratic 
programming, by the interior point method on Lorentz cones, see [AGIIART] . The 
corresponding pocket of programs can be found in http: / /www.mosek.comj . 
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The second difference with respect to Aigorithm (R) is concerned with the 
stopping criterion. In particuiar ([9]) has to be replaced by the following condition 
(|15p . For a given pair (v, A) we have to check in fact the condition 

|Re(i;*, Av)\ < 1 , j = l,...,p. (15) 

4 Computing of the joint spectral radius: the case of nonnegative 
matrices (P) 

Algorithm (R,) can be modified for families of nonnegative matrices to improve 
significantly its efhciency. The corresponding Algorithm (P) has a very similar 
structure, but differs from Algorithm (R) in several key points. Before describing 
the algorithm we need to establish several auxiliary results on operators with an 
invariant cone. 



4.1 Operators with invariant cones. Monotone extremal norms 

Let K he a. convex closed pointed nondegenerate cone with the apex at the origin. 
In the sequel we write K* for the dual cone: K* = {it £ [ mfxeK{u, x) > }. 
According to the Perron-Frobenius theorem, every operator B that leaves a cone K 
invariant has a positive leading eigenvalue Amax = p(S) and K contains a leading 
eigenvector corresponding to this eigenvalue. Any leading eigenvector of B that 
belongs to the cone K will be referred as Perron-Frobenius eigenvector. 

If all operators of the family M share a common invariant cone, then Theo- 
rem [1] on the existence of invariant norms can be slightly sharpened. First, the 
irreducibility condition can be relaxed; second, an invariant norm can always be 
chosen to be monotone with respect to the invariant cone. A function g is mono- 
tone on a cone K if g{x) > g{y), whenever (x — y) £ K . If g is a monotone norm 
defined on the cone K , then it is extended onto M*^ in a standard way: the unit 
ball of that norm is 

{a; e R'^ I ||a;|| < 1 } = cOs {x e K , g{x) < l} . (16) 

All extreme points of the ball defined by p6p are in the cones K and —K. Since 
the norm of any operator A is attained at an extreme point of the unit ball, we 
see that if A leaves K invariant, it attains its norm in the cone K. Thus, 

11^11= max g(/la;). 

xeK, g{x)<l 

In particular, if g is an extremal norm for a family A4, i.e., maxi=i_...^m, \\Ai\\ = p, 
then its extension defined by (|16p is extremal as well. Thus, for families with a 
common invariant cone it suffices to construct an extremal monotone norm g on 
that cone. 

We are going to show that there exists not only extremal, but invariant mono- 
tone norm on K. Recall, that a norm in K is invariant for A4 if 

max{||Aia;||, . . . , ||A^a;||} = p ||a;|| , x e K. 
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To formulate the main result we need some further notation. A hyperplane L C M*^ 
is called a plane of support of a cone K ii LDK 7^ {0} and LOint K = 0. A face 
of a cone is its intersection with some plane of support. For example, a spherical 
cone has only one-dimensional faces (rays); the faces of the cone K = M.'^ are 
coordinate planes: Fi-^...i,. = {a; G IR+ \ Xi^ = ■ ■ ■ Xi^ = 0} , r = 1, . . . , d — 1. A face 
F of a cone K is invariant for an operator A if AF C F. 

Theorem 3 // operators of a family A4 = {Ai, . . . , Am} share an invariant cone 
K and do not have common invariant faces of that cone, then A4 possesses a 
monotone invariant norm on K . 

In the proof of Theorem [3] we use the following lemma from [P3, section 4]: 

Lemma 4 IPS] For any cone K and for any norm on this cone there is a homo- 
geneous continuous function 7(3;) positive on the interior of K such that for every 
operator B leaving the cone invariant we have ||-Ba;|| > 7(a;)||_B||. 

Proof . The case p = is impossible, because in this case the operators must 
have a common invariant face. This fact is simple, and we omit its proof. If p > 0, 
then after normalization it can be assumed that p = 1. Let us first get an extremal 
norm for J\A. Take any e* 6 intisT* and for each n > 1 consider the function 

gn{x) = sup max (e*, Ax\ , x & K . (17) 
Note the following properties of these functions. 

1. For very n we have grL+i{x) < gnix), so the sequence {gn}neN is monotone. 

2. We have 

sup giix) < 00, 

hence, by the monotonicity of the sequence {gn}, all gn are uniformly bounded 
on the unit sphere. To prove that gi is bounded observe that the set L = 
[x G K \ gi{x) < +00 } is either {0}, or K, or a common invariant face of K. 
The latter contradicts the assumption. Assume L = {0}. Then consider the 
compact set S = {x £ K \ (e*, a;) = 1 }. For each j > 1 we define the set 

Vj = < X e S \ max {e*,Ax) > 2 |. 

Thus Vj consists of vectors for which some product A of length j increases 
the value {e*,x) more than twice. If L = {0}, then Llj>iVj = S, and, since 
all Vj are open in S, from the compactness of S it follows that UjLiVj = 
S for some N. This means that for every x £ K there is a product A of 
length at most A'^ increasing the value {e*,x) at least twice. Applying this 
argument successively k times, we conclude that for every x £ K there is 
a product 77^ of length 1^ < kN such that (e*, Ukx) > 2*^ (e*, x), and hence 
ll^Z'fell > 2'° (e*,x) ||e*||~^ (the norm of 77fc is the Euclidean) . Taking the power 
1/lk and the limit as fc — > cx), we obtain p > 2^^^, which is a contradiction. 
Thus, L 7^ {0}, and hence L = K. Thus, each function gn is bounded on S. 
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3. For each n the function Qn is homogeneous, positive (as a supremum of positive 
values) and convex on K, as a pointwise supremum of linear functionals. Thus, 
gn(-) is a norm on K. 

4. For every x ^ K we have gn{Aix) < gn{x) , Ai £ A4. Hence for each n the 
norm gn is extremal. Thus, we have established the existence of a monotone 
sequence of extremal norms. 

5. For every operator A leaving K invariant we denote by 

||A||e* = sup {e*,Ax) 
xes 

the operator norm corresponding to the norm ||a;||e* = {e*,x). Since p = 1, it 
follows that 

max lUlle* > 1. 

Now using Lemmadl we obtain max ||ylx||e* > 7(a;)||a;||e*. This holds for 

Ae^4'' 

each fc, therefore, gnix) > 7(a;)||a;||e* for every x £ K. 

6. Since the sequence {gn{x)}n^N is non-increasing and bounded below, it con- 
verges to some limit function g{x). For every x £ mtK we have g{x) > 
7(a;)||a;||e* > 0. Thus, the function g is convex, positively homogeneous, and 
invariant, i.e., possesses the property 

g{x) = max g{Ajx) . 

J — 1 , .... m 

It remains to show that g is positive on K, in such case it constitutes an invariant 
norm. Since g{x) > for x G intK, we see that the set L = {x £ K \ g{x) = 0} 
lies on the boundary of K. This set is obviously convex, hence it is contained 
on a face of K. Let V be the minimal (by inclusion) face containing L. Since 
AjL C L , j = 1, . . . m, it follows that AjV C V , j = 1, . . .m, which contradicts 
the assumption. Thus, g is an invariant norm, which completes the proof. 

□ 

Now we focus on the case of nonnegative operators, i.e., operators defined 
by nonnegative matrices (which means, with nonnegative entries). Each family 
of nonnegative operators share an invariant cone R^f.. A family of nonnegative 
operators is called positively-irreducible if they do not have common invariant 
faces among the coordinate planes. Applying Theorem [3] to the case K = M^J., we 
obtain: 

Corollary 1 A family of nonnegative positively-irreducible operators possesses a 
monotone invariant norm on M."^. 

Remark 6 The assumption of Corollary [1] is not restrictive, because the general 
case of nonnegative matrices is reduced to the case of matrices without invariant 
coordinate planes. If the matrices possess common invariant planes, then after 
a suitable permutation of the basis vectors all the matrices get a block upper- 
triangular form. The joint spectral radius of the matrices equals to the largest 
joint spectral radius of the blocks. Thus, the problem of JSR computation comes 
to several similar problems with nonnegative matrices of smaller dimensions. A 
fast polynomial procedure to realize this reduction can be found in [JPBll section 
2]. 
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Before describing Algorithm (P) we formulate an analogue of Lemma [T] for 
nonnegative operators. 

Lemma 5 Let a nonnegative operator B have a unique simple Perron-Frobenius 
eigenvalue X with an eigenvector v; let also v* he the Perron-Frobenius eigenvector 
of B* such that (v* , v) = 1. If for some nonnegative operator C one has {v* , Cv) > 
1, then for sufficiently large r the operator B^C has a unique simple Perron- 
Frobenius eigenvalue bigger than A. 

The proof is literally the same as for Lemma [1] □ 
If a cone K C K*^ is fixed, then for a given set Q C M'* we denote 

co_ (Q) = ( co(Q) - K) n K = { X e K \ X = y - z , y e CO (Q) , z e K} . 

If the cone K is not specified, we always assume K = M."^. 

Everywhere below in this section the family A4 is assumed to be positively 
irreducible (see Remark [6]). 



4.2 Algorithm (P) versus Algorithm (R) 

Algorithm (P) ia similar to Algorithm (R) but has some peculiar differences which 
we remark in the sequel. The second has a major computational importance. 

(i) By the Perron-Frobenius theorem the candidate s.m.p. 77 has a nonnega- 
tive leading eigenvalue Amax = p{n), and the corresponding eigenvector be- 
longs to R!j.. We assume that Amax > 0. Hence the main assumption for Algo- 
rithm (R) holds true automatically. 

(ii) The LP problem performed in the loop at step k should be replaced by the 
following: 

max to 
subject to to z < txX 

and X] < 1' > Vx e Vfc 

where we recall that z = Av. 

(iii) In contrast to the case (R), we work now only with nonnegative vectors, and 
do not include the vectors —Vi to the polytope Pk- We do not construct a sym- 
metric polytope cOs (V^), but a positive polytope co_ (Vfc). So, our norm will 
have a unit ball co_ (Vfc). This explains the differences between LP-problems 
(O and (fT8l). 

(iv) Condition ([5} in the Stopping criterion should be replaced by the following: 



[v* , Av) < 1 , j = l,...,n 



(19) 
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4.3 Explanations and proofs 

The theoretical base of Algorithm (P) is actually the same as for Algorithm (R). 
Let us only stress the distinctions. First of all, we construct a monotone extremal 
norm on the positive orthant M^j., so the polytopes Pk are in the orthant, and 
are not centrally-symmetric. That is why we do not need to symmetrize the con- 
vex hull of Vfe. Consequently, LP-problem (|18p . in contrast to LP-problem 
does not have extra variables Qx ,x 6 Vk- The other difference is that the equal- 
ity constraint toAv = YlxeVk ^^'"^ becomes an inequality. This means that the 
polytope Pk is not a convex hull of Vfc, but co_(Vfc). Thus there are two ad- 
vantages of Algorithm (P): 1) the number of variables and the number of con- 
straints in the LP-problem is a half of the LP-problem in Algorithm (R) ; 2) the 
polytope CO- (Vfc) is larger than co(Vfc), therefore this algorithm sorts out more 
vertices ("dead branches") at each step, which leads to a lower complexity. In 
practice Algorithm (P) works much faster than Algorithm (R) (see Section[7]and 
Section [9|). 

In the worst case, if the algorithm has not terminated, one gets an approximate 
value of JSR from inequality 

[p(77)]i/" < p{M) < (t^)-i[p(77)]i/". (20) 

Proposition 3 For an arbitrary positively-irreducible family A4 we have t]\j ^ 1 
as /—> oo and N ^ oo in estimate i20\) . 

The proof is the same as for Proposition [1] with the use of Theorem [3] instead of 
Theorem [TJ The proof of the efHciency of the stopping criterion (Proposition d]) is 
also the same as for Proposition [2l 

Proposition 4 Assume the Perron-Frobenius eigenvalue of 11 is unique and sim- 
ple. If the assumption of the algorithm is wrong (i.e., U is not an s.m.p.) then for 
every j = l,...,n condition II19\) is violated at some step. Conversely, if condi- 
tion it jffj) is violated at some step for some j , then 11 is not an s.m.p. 

Remark 7 In Algorithm (P) the family A4 is assumed to be positively irreducible. 
Actually, this was done for the sake of simplicity. The algorithm can be applied to 
arbitrary nonnegative families. If the algorithm terminates after fcth iteration, and 
the set Vfc does not lie in a coordinate plane of a smaller dimension, then Pfc_i 
is an extremal polytope, and 77 is an s.m.p. That condition means that for each 
i = 1, . . . ,d there is a vector from Vfc with strictly positive ith coordinate. This 
simple condition allows us to apply Algorithm (R) to arbitrary family, without 
preliminary checking its positive irreducibility. Nevertheless, if the family A4 is 
reducible, then it is always advisable to factorize it before starting the algorithm, 
because this significantly reduces the dimension (see Remark [6|). Especially as the 
factorization is realized by a fast polynomial routine [JPBl] . 

5 The criterion for finite termination of Algorithms (R), (C) and (P). 

Algorithms (R), (C), and (P) compute JSR by step-by-step constructing a poly- 
tope norm. If the algorithm terminates within finite time, then it produces an 
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extremal polytope, and, hence, proves that the chosen product 77 is an s.m.p. If it 
does not terminate, then it gives upper and lower bounds for JSR that converge 
to the exact value. An important issue is the following: what are the conditions 
for the family A4 and for the product 77, under which the algorithm terminates 
and produces the extremal polytope ? Conditions guaranteeing the convergence of 
the algorithm to an extremal polytope norm have been discussed in jGWZl[GZ2| . 
We give here a further result which is related to those obtained in the mentioned 
papers. 

Certainly, the product 77 must be spectral maximizing for that. This condi- 
tion, however, does not guarantee the convergence of the algorithm. It appears 
that a bit stronger condition solves the problem completely: it is both sufhcient 
and necessary. The product 77 has to be not just maximizing but dominant. To 
formulate the criterion we need some further notation. 

Let A4 = {Ai, . . . , Am} be a given family of operators, 77 = Ad„ ■ ■ ■ Ad-^ be 
some product, which isjiot a power of a shorter product, n > 1. We denote 
A, = [p(77)]-i/"A,, M = {Ai,...,Am} and 77 = Ad„---Ad,. Clearly, the 
spectral radius of any power of 77 or of any power of its cyclic permutation is 1. 

Definition 3 A product 77 € A4^ is called dominant for the family A4 if there is 
q < 1 such that the spectral radius of every product of operators of the normalized 
family M , that is not a power of 11 nor a power of its cyclic permutations, is 
smaller than q. 

Obviously, the dominant product along with all its cyclic permutations are all 
s.m.p., but vice versa. 

The following theorem gives a sharp criterion on the family A4 ensuring that 
our algorithm produces an extremal polytope. 

Tiieorem 4 For each of Algorithms (R), (C) and (P) the following holds: 

the algorithm terminates within finitely many iterations if and only if 11 is 
dominant for A4, and its leading eigenvalue is unique and simple. 

The proof is in Appendix. 

Corollary 2 If a family A4 possesses a dominant product, whose leading eigen- 
value is unique and simple, then it has an extremal polytope. 

Remark 8 Theorem [3] remains true even if we do not apply the stopping criterion 
in the algorithms. 

Remark 9 Theorem [3] implies that if the algorithm terminates within finite time, 
then the family Ad possesses a dominant product. Actually the algorithm ensures 
that a chosen product 77 is dominant. In numerical examples from applications 
(Sections [S]) and from randomly generated matrices (Section [2| most of matrix 
families possess dominant products. 

The assumption on a dominant product allows to exclude a limit spectrum 
maximizing product, that is a matrix in the closure of the multiplicative semigroup 
of M with spectral radus equal to 1 (see [GZ4 ). 

Consider in fact the following example. Let Ad = {A\, A2): 



Ai={l\) and A. = ^(;j) 
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We can prove that p{M.) = 1+ and there is a unique finite spectrum maximizing 

product n = A1A2 such that p{P) = 1 (apart from its cychc permutation A2A1 
and their powers). Nevertheless the product is not dominant. In fact the sequence 

of products Qk = Ai (A1A2) is convergent and such that 



This implies that P is not dominant; the matrix Qoo is indeed a limit spectrum 
maximizing product of the normalized family M. 

Indeeed the algorithms (R) and (P) do not converge when applied to this 
example. However, if we modify the algorithms and remove a vector when it lies 
on the boundary of the polytope Pk-i, then we obtain a finite convergence also 
in this case. 



6 Computing the lower spectral radius. Algorithm (L). 

In this section we describe a method for the exact computation of the lower spectral 
radius of a finite family of matrices. 

6.1 Antinorms on convex cones 

To extend our approach to computing the lower spectral radius, first of all we 
need the notion of extremal norm for this case. One can define it by the inequality 
miiiAieM > p||a;|| , x G K*^. However, simple examples show that such a 

norm may not exist even for very "good" families M (for instance, irreducible fam- 
ilies of positive matrices). The reason is that the function x 1-^ min^jeM 
may not be convex, in which case it is not a norm (in contrast to the situation 
with JSR, when the function x H- max^.g^vi is always a norm). One of 

the ways to generalize the notion of extremal norm for the lower spectral radius 
is to consider concave positive homogeneous functionals on K"* instead of convex 
ones (i.e., instead of norms). However, such functionals do not exist. Indeed, if / 
is concave, then f{x) + f{—x) < /(O) = (homogeneity failure), hence / cannot 
be positive. Nevertheless, extremal concave "norms" can be defined, provided all 
operators of the family A4 share an invariant cone. In particular, this can be done 
for families of nonnegative matrices. As in the previous section, is a convex 
closed pointed nondegenerate cone with an apex at the origin. 

Definition 4 An antinorm is a continuous nonnegative nontrivial (not identical 
zero) concave positively-homogeneous function on a cone K. 

From the concavity it easily follows that an antinorm can vanish only on the 
boundary of K. An antinorm / is called positive if f{x) > for all x G \ {0}. 
So, / is positive, whenever it is positive on the boundary. This is well known 
that a concave function is continuous at each interior point of its domain. Hence 




which is such that p(Qoo) = 1- 
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the continuity condition for antinorms can be relaxed to the continuity on the 
boundary. 

Let us observe some basic properties of antinorms. First of all, every antinorm 
is asymptotically bounded above by every norm in M'^: 

Lemma 6 For any antinorm f and for any norm \\ ■ \\ there is a constant C such 
that 

fix) < C\\x\\, xeK. 

Proof . Since / is continuous, the value C = sup f{x) is finite. Now by the 

xg_ff,||ri;|| = l 

homogeneity the lemma follows. 

□ 

Consider now a family of operators A4 = {Ai, . . . , Am} that share an invariant 
cone K . 

Proposition 5 If for some antinorm f and for a constant A we have f{Aix) > 
A f{x) , X e K , Ai e M, then p > X. 

Proof. Applying Proposition [5] for an arbitrary point e 6 vaXK , \\e\\ = 1, we get 
\\Aa,---Aa,e\\ > C" • • • ^d,e) > A'^ /(e). Thus, 

min \\Aa,---AaA > C-'f{e)\\ 

dk,...,di 

Taking the power 1/fc and the limit as fc — t- oo, we conclude the proof. 

□ 

Definition 5 An antinorm is called extremal if f{Aix) > pf{x), x & K , Ai £ 
M. 

Similar to monotone norms, an antinorm is called monotone if f{x) > f{y), when- 
ever X — y £ K . 
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Theorem 5 For every family of matrices with a common invariant cone K there 
exists a monotone extremal antinorm on K , 

Proof. If p = 0, then any antinorm suffices, ff p > 0, tlien after normalization it 
can be assumed that p = f . Take any e* G int A'* and consider the function 

fix) = inf min {e*,Ax), xeK. (21) 

This function is concave, homogeneous and monotone, and f(^AiX~) > f{x) for each 
Ai £ A4. It remains to show that / is not an identical zero. Consider the compact 
set S = {x G K \ (e*, x) = 1} and for each k > I define the set 

Uk = I X £ S \ min (e* , Ax) < - |. 

L I AeM'' 2 ) 

If / = 0, then Ufc>i = S, and, since all Uk are open in S, from the compactness 
of S it follows that IJfcLi Uk = S for some A''. This means that for every x G 
K there is a product A of length at most N reducing the value (e* , x) at least 
twice. Applying this argument k times, we obtain that for every x £ K there 
is a product Uk of length Ik < kN such that (e*,77fca;) < 2~'^(e*,a;). Note 
that if j; G int A, then there is a constant C that depends on e* and on x, and 
such that ||S|| < C{e*,Bx) for any operator B that leaves K invariant (see, for 
instance 'P^). Thus, ||77fc|| < C2~^ {e*,x) for each k. Taking the power l/h and 
the limit as fe — oo, we get p < 2~^^^ , which is a contradiction. 

□ 

Applying this theorem for the case K = M[J., we obtain: 

Corollary 3 For an arbitrary family of nonnegative matrices there is a monotone 
extremal antinorm on M."^. 

In the algorithm we need the following analogue of Lemma [5] whose proof is 
the same. 

Lemma 7 Under the assumptions of Lemma\5i if for some nonnegative operator 
C one has (v* ,Cv) < 1 then for sufficiently large n the operator B"C has a 
unique simple Perron- Frobenius eigenvalue smaller than A. 

Now we are ready to describe Algorithm (L) of LSR computation for non- 
negative matrices. We use the notation co-i-(A') = co(A') + K = \^x + h \ x E 
CO {X) h e K}, where X is a subset of R'* and if C M'* is a cone. If X if finite, 
then co+(X) will be referred as an infinite polytope. Thus, an infinite polytope is 
a set P + K, were if is a cone and P is a polytope. In the algorithm we always 
assume if = . 

A product il G is called the spectral lowest product (s.l.p.) if [p(i7)]^/" = 
p{Ai). We shall also call it a spectrum minimizing product, but always use the ab- 
breviation s.l.p. to avoid confusion with the spectral maximizing product (s.m.p.). 
By inequality ^ a product il is an s.l.p. iff [p(i7)]^^" < p{M). Proposition [5] im- 
plies that if there is an antinorm /: R+ R+ such that f{Ajx) > [p(i7)]^/"/(x), 
X G R+, then [p(ii)]^''" = p{A4), and / is extremal. The main idea of the algo- 
rithm is to select a candidate ii for s.l.p. (by a reasonable exhaustion) and then 
to prove that it is actually an s.l.p. The proof is by step-by-step constructing an 
extremal infinite polytope, which generates an extremal antinorm. 
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6.2 Algorithm (L) 

Initialization. We have an arbitrary family A4 = {Ai, . . . , Am}- For some (as 
large as possible) I we look over all products 77 of length < I and take one with 
the smallest value [p(77)]^''", where n is the length of the product. We take the 
shortest product possessing this property and denote it as 77 = Ad„ ■ ■ ■ Ad-^ . If 
p(77) = 0, then p{M) = 0, and the algorithm terminates. So, we assume p(77) > 
0. By the Perron- Frobenius theorem 77 has a positive leading eigenvalue, and 
the corresponding eigenvector belongs to R!j.. We normalize the family A4 as 
Ai = [p(77i)] "^/"Ai , M = {Ai}YLi- The leading eigenvalue of the operator 
77 = Ad^ ■ ■ ■ Ad-^ equals to 1. 

Let 77i = 77 , 77i = Adi_-^ ■ ■ ■ Ad-^Ad,-^ ■ ■ ■ Ad^ be a cyclic permutation of 77i, 
i = 2, . . . ,n. We take a Perron-Frobenius eigenvector vi of 77i (if it is not unique, 
take any of them), and 

Vi = • • • Ad^Vi . 

Thus, Vi is a Perron-Frobenius eigenvector of 77i with the eigenvalue 1. 

In case the leading eigenvalue A max = 1 is unique and simple, we also need a 
dual system of vectors: vl the leading eigenvector of 77i normalized as {vl, vi) = 1 
(Remark [1]), and 

= AX ■ ■ ■ A*d._,vl , i = 2,...,n. 

Thus, Vi and v* are the Perron-Frobenius eigenvectors of 77i and 77* respectively, 
and (w* ,Vi) = 1. 

Set fc = 0. We set Vo = Uq = and TZq = \^{vi,Ap)\ i = 

l,...,n; p= l,...,m, p^d^}. 

Main loop 

For fc > 1. We have finite sets Vk-i C K'^ , Uk-i C Vk-i, and TZk-i C 
Uk-i X M- Put Vfe = Vk-i and Uk = 0. 

We successively take all pairs {v,A) G TZk-i- If for a given pair Av = 0, then 
we stop the algorithm, it is inapplicable for this case. If z = Av ^ 0, then we 
solve the following LP problem with variables to and {tx}xeVk'- 

min to 

^ subject to to ^ tx X /qo) 

and E > 1. tx>0 Vxe Vfc. 
xeVk 

The value of the problem, i.e., min to will be denoted by t^^ Thus, for every 

pair {v,A) G TZk-i we have a nonnegative number t^^ which may take value 
+CXD, when the system of inequality constraints has no solution. 

If t^^ < 1, then leave the sets Vk and Uk as they are, take the next 

pair {v,A) G TZk-i and consider problem (|22p for it. 

Otherwise If 1 < Ti, then 
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If the leading eigenvalue of 77 is unique and simple, we apply the following 
stopping criterion: 

Stopping criterion 

We check the condition 

{v*,Av) > 1, j = l,...,n. (23) 

If (I23l) is satisfied, then we set Vk = Vk U {Av} ,Uk = Uk U {Av}, take 
the next pair {v, A) G TZk-i and consider problem (|22p for it. 

Otherwise If (|23p is not satisfied, then our assumption is wrong, U is not 
an s.l.p., and p{M) < 1 (Lemma O. We stop the algorithm and go either to the 
Final step, or back to the Initialization. In the latter case we need to find 
another pretender to s.l.p. The first opportunity is to increase I and to look over 
all products of a bigger length. Lemma [7] provides also a different approach. We 
take an index j, for which {v* , Av) < 1. Applying Lemma [3 we conclude that 

there is r such that Amax(77J ^s, • • • As-^) < 1, where As,^ ■ ■ ■ As-^vj = Av. We 
take the new initial product 77 = 77]" As^ ■ ■ ■ Ag-^ and restart the algorithm. 

End If 

Otherwise If the leading eigenvalue of 77 is not unique or multiple, then we 
do not apply the stopping criterion, and set Vfe = VkU { Av} , Uk = {Av}, 
take the next pair {v,A) G TZk-i and consider problem (|22p for it. 

End If 

The kth step is over, when all pairs {v,A) G TZk-i UTe exhausted. 

If Uk = 0, then p(M) = 1, and so p{M) = [p(77)]^/". The extremal infinite 
polytope is Pk-i = co^ (Vfe) , and the product 77 is an s.l.p. for A4. The algorithm 
terminates after the fcth step. 

Otherwise If Uk ^0, then we set 7v!,fe = Uk x M and go to the (fc + l)st step. 

End If 

End For 

Final step. If the algorithm has not terminated, then we stop it after some N 
steps, denote tN = jnax tr and have the following estimate for the 

lower spectral radius: 

(t^)-i[p(77)]i/" < p{M) < [p(7r)]i/". (24) 



End of Algorithm (L) 

6.3 Explanations and proofs 

The algorithm produces a sequence of embedded infinite polytopes Pi C P2 C 
... such that Pj+i = co^ {AiPj, . . . , AmPj} and Pj C Pj+i for every j. If 
the algorithm terminates after the fcth step, then Pk = Pk-i- The fcth step is 
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actually needed only to ensure that the infinite polytope Pk-i is extremal, i.e., 
Aj Pk-i C Pk-i ,j = 1, . . . ,m. Since ^ Pk-i (otherwise at some step we had 
Av = 0, in which case the algorithm would be stopped by the end of the fcth step), 
the antinorm fk-i{x) = sup | tx G Pk-i} is well-defined on M'^. Since 

fk-i{Ajx) > fk-i{x) for every x, we see that p{A4) > 1. On the other hand, 
piM) < [p(i7)]i/" = 1. Thus piM) = 1 and so p(Ow) = [p(77)]^/". Thus, if 
the algorithm terminates within finite time, then the s.l.p. and the exact value of 
LSR are found. In this case Pk-i is an extremal infinite polytope and fk-i is an 
extremal antinorm. 

Although the extremal antinorm is obtained numerically, the results are actu- 
ally exact, because the algorithm removes only those points v, for which the strict 
inequality fi{v) > 1 holds, where fi is the antinorm generated by the current 
polytope Pi = CO+ (Vi). 

If the algorithm does not terminate within finitely many steps, then we have es- 
timate (|24p to get an approximate value of the LSR. The right hand side inequality 
is obvious, the left hand side is equal to min inf fi\[-i{Ajx) < p{A4), 

j = l,...,m /„_i(a:) = 1 

from which the estimate follows. 



6.4 Efficiency results for Algorithm (L) 

Let us start with the stopping criterion. If the stopping criterion is applicable 
(i.e., the leading eigenvalue of 77 is unique and simple), then it always determines, 
whether 77 is an s.l.p. or not. 

Proposition 6 Assume the Perron- Frobenius eigenvalue of U is unique and sim- 
ple. Ij the assumption of the algorithm is wrong (i.e., 77 is not an s.l.p.) then for 
every j = 1, . . . ,n condition \23\) is violated at some step. Conversely, if condi- 
tion S23\) is violated at some step for some j, then 11 is not an s.l.p. 

Proof. The sufficiency follows from Lemma [T] To^ow the necessity, we assume 
that p{M) < 1. Then there is a product C G M" such that p{C) < 1. This 
yields — > 0, and hence C^vi — > as r — > cxd. Consequently, for every j one has 
{v* , C^vi) < 1, whenever r is large enough. Since the point C^vi belongs to the 
infinite polytope Prs, we see that inf (v*,C^vi) < 1. This infimum is attained 

at some vertex v of Prs , hence that vertex violates condition (|23p . 

□ 

Now let us analyze estimate (|24p . In contrast to the algorithms for JSR com- 
putation, the lower bound {tis[)~^[p{n)]^^" may not converge to p at all, even if 
the family M is positively irreducible. There are simple examples already in the 
dimension (7 = 2. The reason is that some product A G A^' may have the leading 
eigenvector v on the boundary of the invariant cone i.e., have some zero en- 
tries. If the corresponding leading eigenvalue is unique and simple, then for each 
j = 1, . . . ,n the sequence A^vj converges to tjW as r —)■ oo, where tj > depends 
only on j. This means that some vertices of Pk approach closer and closer to the 
boundary as A; oo. In this case Algorithm (L) is useless: it gives neither an 
extremal infinite polytope nor good lower bound for p. We suggest two methods 
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to avoid this situation. The first one is to impose a second invariant cone assump- 
tion. Second, is to enlarge the invariant cone by adding extra directions. Let 
us begin with the first one. 



6.5 Convergence of the algorithm. Case 1; the second invariant cone. 

A cone K is embedded to a cone K if {K \ {0}) C int K. The pair {K, K) will be 
refereed as an embedded pair. 

Definition 6 An embedded pair is invariant for a family M. if both its cones are 
invariant for this family. 

An embedded invariant pair of cones for a given family will be called an invariant 
pair. For a family A4 of nonnegative operators (i.e., operators defined by nonneg- 
ative matrices) we say that K is a second, invariant cone if {K.W^) is an invariant 
pair. So, the cone K is embedded in W^. and invariant for A^. A simple sufficient 
condition for the existence of a second invariant cone is the so-called eventual 
positivity of a matrix family. 

Definition 7 A nonnegative family A4 is called eventually positive if there is k 
such that all matrices of the family A4^ are positive for all r > k. 

In particular, if all matrices of M. are positive, then M. is eventually positive. The 
following trivial fact clarifies the notion of eventual positivity: 

Lemma 8 A family M. is eventually positive iff its matrices have neither zero 
columns nor zero rows, and there is k such that all matrices of M'^ are positive. 

Lemma 9 Every eventually positive family possesses a second invariant cone. 

Proof . A conic hull of the set U^eM''-^^^ where K = is the second invariant 
cone for M. 

□ 

The key property of embedded pairs is formulated in the following lemma. 

Lemma 10 If {K, K) is an embedded pair, then every antinorm on K is contin- 
uous and strictly positive on K. 

Proof . Any nonnegative concave function, which is not an identical zero, is con- 
tinuous and positive at any internal point of its domain. 

□ 

Now we can prove the existence of invariant antinorms in the interior cone. 

Definition 8 Let a family A4 has a common invariant cone K . An antinorm 
f : K ^ K+ is called invariant, if there is a constant A > such that 

min f{Ajx) = Xf{x) , x £ K. 

j=l,...,m 

Theorem 6 // the family Ai possesses an embedded pair {K,K), then it has a 
positive monotone invariant antinorm on K. For any invariant antinorm on K 
we have A = p{Ai). 
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Let us recall that for extremal antinorms we have f{Ajx) > px. It 

becomes an equality if the antinorm is invariant. 

By Theorem [S] an extremal antinorm always exists, whenever the operators 
share an invariant cone. For the invariant antinorm this is not the case. There 
are simple examples of irreducible pairs of nonnegative 2 x 2-matrices that do 
not have an invariant antinorm. So, the embedded pair assumption is essential in 
Theorem (6] In the proof of Theorem [6] we use the following simple fact from |P31 
section 4] . 

Lemma 11 For every pair of embedded cones {K,K) and for every norm in R*^ 
there is a constant 7 such that for any operator B with these invariant cones and 
for each x (z K one has ||-Ba;|| > 7 ||i3|| ||a;|| . 

Proof of Theorem [6]. Let / be an invariant antinorm. For each k we have 
miuyig^fc f{Ax) = \^ f(x). On the other hand, combining Lemmas [6l and [TT1 we 
see that for any operator B that preserves the cones K and K one has 

C4B\\\\x\\ > fiBx) > Co'\\B\\\\x\\, xeK, 

where Co does neither depend on B nor on x. 

A nontrivial concave nonnegative function is strictly positive on the interior of 
its domain. Therefore, /(x) > for all re G K\{0}. By the compactness argument 
it follows that f{z) > Ci||2:|| for any z G K, where Ci > does not depend on z. 
Hence, applying Lemma[TT] we conclude that f{Bx) > d||-Ba;|| > d7||i3|| 

Therefore min^^g^i- \\Ax\\ x A*^, and therefore p = X. Now let us prove the ex- 
istence of an invariant antinorm. It sufhces to consider the case p{A4) = 1. By The- 
orem [5] there is a monotone extremal antinorm /o, for which mmAeM fo{Ax) > 
fo{x) , X 6 K. Let fj{x) = minAeM fj-i{Ax) , j 6 N. This is a nondecreasing 
sequence of antinorms. If for some x £ K we have fj{x) -\-oo as j — > 00, 
then this holds for all nonzero x £ K, and hence p{A4) > 1. Thus, the sequence 
{fj}jeN is bounded, hence it converges pointwise to some function /, which is a 
monotone invariant antinorm. 

□ 

Applying now invariant antinorms we can prove the convergence results for 
Algorithm (L). We start with inequality (|24p . 

Proposition 7 // the family A4 possesses a second invariant cone A' C M+, and 
vi G K , then for estimate ^24\ ) we have tjv ~^ 1 ffls Z — > 00 and N — > 00. 

Thus, if A4 has a second invariant cone, then Algorithm (L) is always applicable 
for, at least, approximate computation of LSR. It either finds the value of LSR or 
provides lower and upper bounds for it; those bounds are arbitrarily close to each 
other, whenever both / and N are large enough. 

In the proof we use Dini's theorem (see [Ru', theorem 7.13]) and the following 
analogue of the Minkowski norm for concave functionals. We call a convex closed 
set D C 1R+ admissible if it does not contain the origin, and if with every point x G 
D it contains all points y > x. In particular, all infinite polytopes not containing 
the origin are admissible. The Minkowski antinorm associated to an admissible set 
D C ]R+ is defined as foix) = sup{t"^ \ t > 0, tx e D} , x e M+. For any 
admissible set D the function /u is a positive monotone antinorm on W^. 
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Proof of Proposition IT] Assume first that p{A4) = 1, i.e., that 77 is an s.l.p. 

The algorithm produces the infinite polytopes {-PfcjfcgN such that Pk C Pk+i = 
CO+ {AiJ-fc, . . . , AmPk}- All their vertices are in K, because vi G K. By Theorem[6] 
there is a positive invariant antinorm fo on K, for which 

foix) = min fo{Ajx) , x e K. 

j = l,...,m 

Therefore, for all vertices v of the polytopes Pk one has fo{v) > fo{vi). Since fo is 
positive on K, by the compactness argument it follows that fo{x) > C\\x\\ , x ^ K, 
where C > is a constant. Whence, all polytopes Pk are uniformly separated from 
zero: they do not intersect the ball of radius C||?;i|| centered at the origin. Conse- 
quently, the sequence {/fc}fceN of Minkowski antinorms generated by the infinite 
polytopes {Pfc}fcgN is non-decreasing and bounded. So, it converges pointwise to a 
positive monotone antinorm /. By Dini's theorem [Rul theorem 7.13], this conver- 
gence is uniform on the set S = {x ^ K \ f{x) = 1}. Thus, fk(x) — > 1 uniformly 
for a; G 5, as fc — > oo. Hence, there is Ne such that /iv-i(x) > 1 — £ for all x £ S, 
whenever > iVe. Hence, (<jv)"^ = inf^^gp^^ /jv-i(a;) > inf^gs /at-i (a;) > l-£, 
which completes the proof for the case p{M) = 1. The transfer to the general case 
is realized in the same way as in the proof of Proposition [1] 

□ 

Corollary 4 Each of the following conditions is sufficient for the convergence 
t]\f ^ 1 as I ^ oo and N ^ oo: 

1 ) the family A4 is eventually positive; 

2) the family A4 has a second invariant cone, and the leading eigenvalue of 11 
is simple. 

Proof . Observe that if a family is eventually positive, then every product 77 of 
its matrices has a unique simple largest by modulo eigenvalue. This eigenvalue is 
positive, and the corresponding eigenvector is strictly positive. To see this note 
that the matrix 77 is obviously primitive (i.e., it is nonnegative and some power 
77*^ is strictly positive). A primitive matrix always has a unique simple largest by 
modulo eigenvalue, which is positive, and the corresponding eigenvector is strictly 
positive [HJI chapter 8]. If 7V4 is eventually positive, then for each A £ Ad'^ and 
every v G 1R+ the vector Av belongs to the interior cone K, which is a conic hull 
of the set |J^^ Ad^. ■ ■ ■ A^-^ Hence, K contains the leading eigenvector of 

any product 77 of matrices from A4. 

If A4 has a second invariant cone K, then, by the Perron-Frobenius theorem, K 
contains some of the leading eigenvectors of 77. If the leading eigenvalue is simple, 
then vi G K. 

□ 



6.6 The criterion for finite termination of Algorithm (L). 

Now we are ready to prove a sharp criterion ensuring that Algorithm (L) produces 
an extremal infinite polytope. It looks similar to Theorem |3] and use the notion of 
under-dominant product. 
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Definition 9 A product 77 £ is called under- dominant for the famibj^Ad if 
there is p > 1 suchj,hat the spectral radius of every product of operators of M , that 
is not a power of U nor a power of its cyclic permutations, is bigger than p. 

Tiieorem 7 Assume the family A4 is eventually positive. Algorithm (L) termi- 
nates within finitely many iterations if and only if 11 is under- dominant for A4. 

The proof of Theorem [7] is in Appendix 1101 

Remark 10 If the family J\A is eventually positive, then every product 77 G A^" is 
a primitive matrix, i.e., some of its powers is positive. 

This is well known that the leading eigenvalue of a primitive matrix is always 
unique and simple (see, for instance, |HJI chapter 8]). That is wliy in TheoremOwe 
do not need the uniqueness and simplicity of the leading eigenvalue assumption, 
in contrast to Theorem |31 

Remark 11 Let us stress again that Algorithm (L) can produce extremal polytopes 
for nonnegative families that have no second invariant cone or not eventually 
positive. We will see some examples in Section [3 and [8l The only difference is 
that we have not succeeded in finding a reasonable criterion for that case. As for 
Theorem[71 the eventual positivity assumption is essential and cannot be omitted. 

Corollary 5 If an eventually positive family A4 possesses an under- dominant 
product, then it has an extremal infinite polytope. 

6.7 Case 2. Modification of Algorithm (L) 

In the previous subsection we showed that if the family M has a second invariant 
cone (in particular, if this family is eventually positive), then the algorithm is 
always applicable. It may converge within finite time, in which case it produces 
the extremal infinite polytope and finds the exact value of LSR. By Theorem[7]this 
happens precisely when the product 77 is under-dominant. Otherwise, if it does 
not converge, it produces upper and lower bounds in (|24p that both tend to p{M) 
as A'' — > oo (Proposition [71). If M does not have the second invariant cone, then 
the algorithm can be applied as well, but in some cases it may not converge to the 
value of LSR. This happens, for instance, when there is a product A £ JvV , whose 
leading eigenvalue Amax is unique and simple, and the corresponding eigenvector v 
has some zero entries. If [Amax]^^'' > [pin)]^^"^ then, for each point vi produced by 
the algorithm we have A^Vi — >• tiV as s — > cxd, where the sequence ti >Q diverges. 

Whence, some vertices of the polytopes Pk converge to the boundary of as 
fc — > oo. In this case the algorithm does not terminate within finite time, since new 
vertices Vi will always appear (closer and closer to the boundary of M^j.). Moreover, 
the ratio tjsi va (|24|1 may not converge to 1, and the algorithm becomes useless. 
We suggest the following modification of the algorithm for this case, which often 
leads to the precise values of LSR. 

Assume the algorithm has not terminated after A'' steps. Take some small 
(5 > and find all vectors v G Vjv such that < ^, where 

'^^min and i^max are 

respectively the smallest and the largest entry of ti. To any such a vector v we 
associate a vector h such that /i' = —e if -^p — < ^, and /i' = 1 otherwise (we 
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write for the qth entree of the vector x) . The parameter e > is chosen to be 
smaU and the same for all h. The finite set of all vectors h will be denoted as H. 
We also denote by e the vector of ones. 

For every j = 1, . . . ,m and for every h & Ti we solve the following LP problem: 



max £e 

subject to Aj h > tf,e + h (25) 

hen 

and th > 0, Mhen. 

If for some j and h we have < 0, then for the chosen values of 5 and £ 
the modification is impossible. We can try smaller values. If ie > for all j and 
h £ H, then we restart our algorithm with the same product 77 and with the only 
modification: LP problem (|22|1 is replaced by the following LP problem, where 
z = Av: 



min to 

subject to ^0 > Yl tx X + Yl thh 

xeVk hen 

and Yl tx>l, tx>0 yxe Vk, 



Explanation, li te > in LP problem (|25p for every j = 1, . . . ,m and for every 
h £ H, then the cone K-u = {x + Yhen ^''^ I ^ ^ "^1- ' ^h > , h £ H} is 
invariant for the family A4. If the algorithm terminates after the fcth step, then 
the set Pk-i = CO (V^v) + K-h is an extremal infinite polytope for the family M, 
i.e., APk-i C Pk-i for all A G M. Therefore, it defines an extremal antinorm in 
the cone K-u, and hence p{A4) = 1. Thus, we replace the invariant cone M.'^ by a 
wider invariant cone K^, which covers those vertices v G Vjv that come too close 
to the boundary of R+. In many practical cases this trick makes the algorithm 
converge within finitely many steps. We use it in the proof of Theorem |8] in ii8.ll 

We consider in the sequel several numerical examples of implementation of our 
algorithms, and start with the simplest case of nonnegative 2 x 2-matrices. In Ex- 
ample[T]Algorithm (P) finds the JSR of two matrices, in Example[2]Algorithm (L) 
finds the LSR of another pair of matrices. The aim of those examples is to show 
how the algorithms work. Then in Section [8] we apply our algorithms to matrices 
of bigger dimensions (up to d = 50) arising in various problems of combinatorics 
and number theory. Finding the exact values of JSR and LSR we prove, in par- 
ticular, several previously stated conjectures in combinatorics and number theory, 
and disprove one conjecture on Pascal's rhombus. Further, in Section [9l we show 
the statistics how our algorithms work for randomly generated matrices of various 
dimensions. For all randomly generated families the algorithms found the exact 
values of JSR and of LSR (the latter is in the case of nonnegative families) . 



(26) 
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7 Illustrative examples 

We consider liere some simple examples showing the flow of the algorithms we 
have presented both for the joint spectral radius and for the lower spectral radius 
of a nonnegative set of matrices. 

Example 1 .Computation of the joint spectral radius. Consider a family 
M = {Ai,A2}: 

with b = 9/10. Looking over all matrix products up to some length, we make a 
guess that U = Ai A2 is a spectrum maximizing product. We have [p(77)]^^^ = 
Vb^^^. Applying then Algorithm (P) to the family M = [p{n)]~'^^'^M we 
obtain at step zero a unique leading eigenvector vi of the product 77i = 77, 
and V2 = A2V1 — the leading eigenvector of the cyclic permutation 772: vi = 
(1 , ^|=i) and 



V2 



A2V1 = Vb( ^ ^ , 1) = (0.586318522 , 0.948683298) 



(all the values are rounded to the ninth decimal). At the first step we get one new 
point 

V3 = Aivi = -^fl , = (1.054092553, 0.402627528), 

y/b\ y/5 + lJ ' 

and the other point V4 = A2V2 is "dead" , because it belongs to the interior of 
P\ = CO- {vi,V2,V3}, i.e., solving LP problem (|18p for the point v = V2 and for 
A = A2 we get t^^^ > 1. Thus, after the first step Vi = {vi,V2,V3} and 

Ui = {V3}. 

At the second step we solve LP problem (|18|1 for the pairs (ws, Ai) and (w3, A2) 
and find that the values t, ^ , and t, t , are both bigger than 1. This means 

^ {y3,Ai} {V3,A2} ^ 

that the points A1V3 and A2V3 are both internal to Pi. Therefore, AjPi C Pi , j = 
1,2, and so Pi is an extremal polytope (see Figure [2|). 

Thus, the algorithm terminates after the second step, and p(A^) = [p(77)]^^^ = 

The cyclic tree of this algorithm is plotted in Figure O In Figure |3] we plot the 
points {AiVi}^^i (in red) and {A2Ui}f=i (in blue). 

Example 2 .Computation of the lower spectral radius. Let M = {^1,^2} 
with 

/70\ , /24 

We prove that the product 77 = Ai A2 {Ai A2)^ is spectrum minimizing, and 

hence the LSR p{M) equals to p(ny/^ = (4 (213803 + ^44666192953))^^^ = 
6.009313489. . .. 
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1.0 




1.0 



2.0 



Fig. 2 The extremal polytope Pi of Example ^ The starting leading eigenvector vi of 77 is 
indicated in red. 



Fig. 3 The cyclic tree of Example [T] The root {vi, V2} is in red, the alive leaves are in green, 
the dead leaves are blue. 



Starting Algorithm (L) we define at zero step M = [p(77)] ^^^M , 77 = 
[p{n)]~^n and get eight points vi, . . . ,vs starting from the leading eigenvector 



of 77. Thus, V2 = A2V1, V3 = A1V2, V4 = A1V3, V5 = A2V4, ve = A1V5, 
V7 = AiVQ, vs = A2V7. At the first step we solve eight LP problems (|22|) and 
get the only new alive vertex vq = A2VQ. 

The other seven new vertices are "dead leaves": they belong to the interior 
of the infinite polytope 7-*! = co+ {vi}\^i, for each of them the value t^^ of 

problem (|22|1 is smaller than 1. Thus, after the first step we have Vi = {i;i}f=i 
and Hi = {fg}. The next step produces two new vertices A^vq , A2Vg and they are 
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1.0 2.0 

Fig. 4 The extremal polytope of Examplc^and the transformed vectors {AiVi}^_-^ (in blue) 
and {A2Vi}^^-^ (in red). 




Fig. 5 The cyclic tree of Example [2] (the red root, green alive leaves and blue dead leaves). 



both dead. Thus, the algorithm terminates after the second step, and Pi is the 
extremal infinite polytope. 



8 Applications 



We consider four applications of Algorithms (P) and (L) to various problems of 
combinatorics, number theory, and theory of formal languages. Each problem is 
reduced to computing JSR or LSR of some families of nonnegative matrices, which 
we are able to solve exactly in the sense specified in previous sections. 
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.0^ 



1.0 2.0 

Fig. 6 The extremal infinite polytope in Example [2] (the starting eigenvector vi of 77 is red). 




1.0 2.0 

Fig. 7 The extremal infinite polytope Pi and the transformed vectors {AiVi}^_-^ (in red) and 
{^2fi}i=i (in blue) of Example^ 



8.1 The asymptotics of the number of overlap-free words 



The problem of counting of overlap-free binary words was intensively studied in 
the literature (see the recent survey [Be] ). In [C] and then in [JPB2] this problem 
was reduced to computing JSR and LSR of two special nonnegative 20 x 20- 
matrices. Those values were computed approximately, and two conjectures were 
stated about their exact values |JPB2] . Now we prove both those conjectures by 
applying Algorithms (P) and (L). 
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A binary word, i.e., a finite sequence of zeros and ones, is caiied overlap-free 
if it does not contain a subword of tfie form xaxax, where x G {0, 1} and a is a 
word. In 1906 Thue proved that there are infinitely many such words. A natural 
problem, which was analyzed in many papers, is to estimate the total number Un 
of overlap-free words of length n. In 1988 Brlek showed that Un > 3 n — 3, on 
the other hand Restivo and Salemi in 1985 proved the polynomial upper bound 
Un < Cn^, where r = log(15) w 3.906. This result was sharpened successively 
by Kfoury (1988), Kobayashi (1988), and Lepisto (1995) to the value r = 1.37. On 
the other hand, Kobayashi (1988) showed that Un > Cn^'^^^ (see [Be for the 
corresponding references and historical overview). So, the number of overlap-free 
words grows faster than linearly. A natural question arises, whether u„ x n'' 
for some 7 G [1.155, 1.37] . Cassaigne [C] showed that the answer is negative. He 
introduced the lower and the upper exponents of growth: 

Q = sup {r I 3 C> 0,u„ > Cn''}, (27) 
/3 = inf {r I 3 C7 > 0,u„ < Cn''}, 

and proved that a < f3. Moreover, he established that the numbers Un can be 
computed as sums of variables that are obtained by certain linear recurrence re- 
lations. This led to the following bounds: a < 1.276 and (3 > 1.332. The next 
improvement is due to Jungers, Protasov and Blondel |JPB2j . who, showed that 
a = log2 p(Ai, ^2) and /3 = logj p(Ai, A2), where yli and A2 are special 20 x 20- 
matrices with nonnegative integer entries, which are reported in the AppendixllOl 
In [JPB2j the authors introduced new algorithms for estimating the joint and 
lower spectral radii based on the convex programming; by means of these algo- 
rithms they derived the following bounds: 

1.2690 < a < 1.2736 and 1.3322 < ,9 < 1.3326. (28) 

This allowed the authors to make the following conjectures on the precise values: 

Conjecture 1 pPB^ 

The s.m.p. for the family A4 is A1A2, and P = \ log2 p{AiA2). 

Conjecture 2 [JPB2] 

The s.l.p. for the family M. is ^1^2°, and a = ji log2 p[A\A^). 

Algorithms (P) and (L) now make it possible to prove both these conjectures. 

Theorem 8 For the upper and lower exponents of growth of the function Un one 
has: 

a = -L iog2 p[AiAl°) = 1.273553265. . . . 

^ = i log2 p{AiA2) = 1.332240491 .... 

Thus, both the upper and the lower exponents of asymptotic growth of the overlap- 
free words can be found precisely. To prove Theorem [8] it suffices to present the 
corresponding extremal polytopes. Since the matrices A1A2 are nonnegative, one 
can apply Algorithm (P) for the JSR computation. The candidate for s.m.p. is 
77 = AiA2. The algorithm terminates having performed fc = 10 steps. Thus, 
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p{Ai,A2) = y^piMM) = 2.517934040.... The extremal polytope Pg has 54 
vertices. 

To compute the LSR we apply Algorithm (L). The candidate for s.l.p. is 
77 = AiA^'^ . However, performing fc = 10 steps we see that the algorithm does not 
converge. There are two sequences of vertices Vi G Vk that approach to the bound- 
ary of the positive orthant i.e., those points have very small entries on some 
positions. Therefore, we apply the modified version of the algorithm (!j677|). Taking 
5 = 1/200, we see that one sequence have entries of index q £ h = {5, 10, 17, 18} 
smaller than S, the other sequence have entries of index q £ I2 = {7,8,15,20} 
smaller than 6. We take e = 1/4 and 7i = {/ii,/i2}, where hj G K^°, the gth 
entry of hj is —s if g G Ij and 5=1 otherwise, j = 1, 2. Solving LP problem (|25p 
for aW h £ Ti , A £ {Ai, A2} we obtain te > for each of those four problems, 
and hence = {a: + tihi + ^2^2 | x G M^j. , ti,t2 > } is a common in- 

variant cone for Ai,A2. Now we apply again Algorithm (L) with the cone Ky^ 
instead of M^j., i.e., replacing LP problem (|22|l by (|26p . The modified algorithm 
converges: it terminates after fc = 15 steps. Thus, 77 = 711^2'^ is an s.l.p., and 
p{Ai,A2) = [p(^iA2°)]^^" = 2.417562630. ... The extremal infinite polytope 
7-'i4 has 104 vertices. 

In order to give the formal proof of the theorem it is sufficient to simply provide 
the list of vertices of the corresponding extremal polytope Pg (to compute JSR) 
and of the extremal infinite polytope 7-'i4 (to compute LSR). They are given in 
the Appendix [TU] 



8.2 The density of ones in the Pascal rhombus 



The Pascal rhombus is related to the Pascal triangle, with the only difference 
that each element equals to the sum of four previous elements rather than two 
(see |GKMT| for definitions and basic properties). The elements of the Pascal 
rhombus arise from linear recurrence relations on polynomials. The sequence of 



polynomials {pn} is defined as po{x) = 1 , pi{x) 



+ X + 1 and pn{x) 



{x + X + l)pn-i{x) + X pn-2{x) ,n > 2 . This leads to a recurrence relation for 
the number Wn of odd coefficients of p„ |FSBj . The asymptotic growth of Wn as 
n — > cxD is characterized as follows: 



lim sup 



log Wn 

log n 



log2 P 



log w„ 



lim inf 

n-KX3 log n 



log2 P, 



where p and p are respectively the JSR and LSR of matrices Ai, A2 defined as 



Ai 



/O 1 0\ 
1 2 

10 1 

\0 2 1/ 



A2 



/I 2 0\ 
2 1 
110 


yo 1 0/ 



(29) 



This is shown easily that p = 2, and the main difficulty is to compute p. 
Conjecture 3 ,!^/ For the set of matrices ^ we have p = = 1.61803. 
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This conjecture is quite natural, because the golden number appears in many prob- 
lems of combinatorics. An approximate computation of p{Ai,A2) given in |PJB| 
provided the following estimate: 

1.6180 < p < 1.6376, (30) 

which rather confirms Conjecture [S] The upper bound in (|30p is obtained by the 
product AlAl, for which = 1.6376.... 

Applying Algorithm (L) we find the precise value of p. It follows that Conjec- 
ture [3] is not true, and the product Ai A2 is actually an s.l.p. 

Theorem 9 For the family M = {Ai,A2} one has p{M) = p^^'^{AlA2) = 
1.6376.... 

This solves the problem of asymptotics of the sequence Wn and disproves the 
golden number conjecture. 

Proof of Theorem [5] To prove the theorem it is convenient to apply Algorithm (L) 
to the transposed family Ai* = rather than to original family (|29p . Of 

course, p{M*) = p{M). We take 77 = {AlfiAl^f as a candidate for s.l.p. The 
leading eigenvector of 77i = 77 is 

VI = (0.39925900 , 0.95496725 , 0.79851800 , 0.90909090 , 1) . 

At step zero we get six vertices: V2 = A2V1, = A2V2, V4 = W5 = 

Alv4, VQ = Alv5. The first step gives two new alive vertices: V7 = Alv3 and 
Us = A2VQ. The algorithm terminates after the second step, since we compute 
that the images of the new vertices V7,vs lie inside the infinite polytope 7-"! = 
co-|-{ui, . . . , vg}. Thus, we obtain the leading eigenvectors of all cyclic permutations 
of 77 plus two additional points (^7 and ws). One can check directly that the infinite 
polytope Pi = co-|-{t;i, . . . ,^8}, where the vertices {vi}i—i are described above, is 
extremal for M* , i.e.. A* Pi c Pi , j = 1,2. 

Thus, the polytope T'l is extremal, and the product {Al)^{A2)^ is an s.l.p. 
In Figure [8] we see the cyclic tree of Pi. Moreover, by Theorem [7] the product 
{Ai)^{A2)^ is under-dominant. Hence, the transpose product /I2 ^1 is under- 
dominant for A4 , and therefore so is the product Ai A2 being its cyclic permuta- 
tion. 

This concludes the proof. 

□ 



8.3 The Euler binary partition function 

For an arbitrary integer r > 2 the Euler binary partition function b{k) = 62 ("r, k) 
is defined on the set of nonnegative integers k as the total number of different 
binary expansions k = J^JLq dj^-' , where the "digits" dj take values from the set 
{0, . . . , r — 1}. The asymptotic behavior of b{k) as fc — > 00 was studied in various 
interpretations by L. Euler, K. Mahler, N.G. de Bruijn, D.E. Knuth, B. Reznick 
and others (see |P3] for the corresponding references). For even r = 2n, as it was 
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1^8 



Fig. 8 The vertices {yi\^^_^ for LSR computation of the Pascal rhombus. 







JSR 






LSR 




r 


# its 


# vertices 


P 


s.m.p. 


# its 


# vertices 


P 


s.l.p. 


7 


5 


8 


3.511547 


Ai 


6 


14 


3.491891 


A1A2 


9 


6 


18 


4.503099 




5 


17 


4.494492 


Ax 


11 


5 


14 


5.505892 


Ai 


7 


24 


5.497042 


A1A2 


13 


5 


16 


6.502167 


Ai 


7 


28 


6.498946 


A1A2 


15 


7 


40 


7.500106 


A1A2 


6 


23 


7.499841 


Ai 


17 


7 


40 


8.500057 




6 


30 


8.499904 


Ai 


19 


7 


24 


9.500423 




8 


46 


9.499789 


AxA2 


21 


6 


28 


10.500373 


Ai 


8 


50 


10.499813 


A1A2 


23 


8 


52 


11.500053 


A1A2 


6 


31 


11.499894 


Ai 


25 


9 


34 


12.500059 


Ai 


8 


58 


12.499971 


A1A2 


27 


8 


60 


13.500030 




7 


37 


13.499938 


Ai 


29 


9 


66 


14.500009 


A1A2 


8 


43 


14.499982 


Ai 


31 


9 


30 


15.500001 


Ai 


10 


34 


15.499999 


A1A2 


33 


11 


36 


16.500001 


A1A2 


10 


55 


16.499999 


Ax 


35 


8 


52 


17.500007 


Ai 


18 


102 


17.499997 


A1A2 


37 


8 


54 


18.500012 


Ai 


10 


113 


18.499994 


A1A2 


39 


10 


112 


19.500003 


A1A2 


8 


59 


19.499994 


Ai 


41 


9 


78 


20.500005 


Ai 


11 


120 


20.499997 


A1A2 



Table 1 Computation of the JSR and of the LSR, for the Euler partition function matrices. 



shown in |Rj , one has h{k) x fc'"^^ " . For odd values of r the asymptotic behavior 
of h(k) is more compUcated and has been studied in m and [P3j . Denote 

pi = liminf log 6(fc)/ log fc; p2 = hmsup log 6(fc)/ log fc . (31) 

In |P3j it was proved that pi = log2p(Ai,A2) and p2 = logj p(Aij42), 
where Ai,A2 are (r — 1) x (r — l)-matrices defined as follows: {As)ij = 1 if 
2 — s<2j — i<r — s + 1, and {As)ij = otherwise (for s = 1, 2). For example, 
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for r = 7 we have the following 6 x 6-matrices: 



/I 1 1 1 0\ 
1110 
11110 
1110 
1111 

\0 1 1 1/ 



A2 = 



/I 1 1 0\ 
11110 
1110 
11110 
1110 

\0 1 1 1 1/ 



In |P3| the following conjecture was made: 

Conjecture 4 For every odd r one of the two products A\ and A1A2 is an s.m.p. 
and the other is an s.l.p. 

The case r = 3 was carried out earlier in the work [R], for r = 5,7,9,11,13 
Conjecture m was proved to hold true in |P3j . Algorithms (P) and (L) make it 
possible to prove this conjecture for many more odd values (in particular we did 
the computation for r < 41). 

The results are listed in Table [TJ The first column is r (where we recall that 
the dimension of matrices Ai,A2 is r — 1), the second column is the number k 
of iterations necessary to Algorithm (P) for terminating, the third column is the 
number of vertices of the extremal polytope Pk-i, the fourth one is the value 
of JSR rounded to the sixth decimal, and the fifth one is the s.m.p. The right 
hand side of the table presents analogous informations for the LSR computation 
by Algorithm (L). 

We see that Algorithms (P) and (L) demonstrate a good efficiency. Even for 
large dimensions of the matrices Ai , A2 the total number of iterations k never 
exceeds 18 and the number of vertices of the extremal polytope Pk-i is at most 
120. Let us remark that the binary matrices A\,A2 of the partition function are 
rather inconvenient for our algorithms, because of a very small gap between JSR 
and LSR. For instance, for r = 33 the distinction between the JSR and LSR is less 
than 0.00002%. Therefore, all products of Ai and A2 of some length k have almost 
the same spectral radii. This is why we would expect Algorithms (P) and (L) to 
need a large number of iterations. On the contrary, they just need 11 and 10 
iterations respectively. 

Actually Algorithms (P) and (L) work also for higher dimensions and Con- 
jecture [Dean be proved for larger r. For instance, if r = 51, then Algorithm (L) 
needs A; = 15 iterations and produces an extremal infinite polytope P14 with 135 
vertices. 

In the next section, as further example, we consider ternary expansions and 
show that also in this case we can compute the significant measures. 



8.4 The Euler ternary partition function 

The LSR and JSR appear in the problem of asymptotics of the Euler partition func- 
tion on the arbitrary base, not only for binary expansions. For instance, the ternary 
partition function b{k) = hz{r,k) is the total number of different ternary expan- 
sions k = X!°^o dj"3-' , where the "digits" dj take values from the set {0, . . . , r — 1}. 
The largest and the smallest exponents of growth of h{k) &s k ^ ozi are defined by 



44 



Nicola Gugliclmi, Vladimir Protasov 




Fig. 9 The vertices {vi}j^-^ for JSR computation of the ternary Euler partition function. 



formulas similar to (|3ip through the LSR and the JSR of three special binary ma- 
trices Ai,A2,A3 (see [P3j for details). In [PJBj the authors analyze the example 
with r = 14, where the matrices of = {Ai, A2, A3} are 

/I 1 1 1 1 0\ 
111110 
11110 
A2 = 111110 
111110 
11110 

Vo 1 1 1 1 1/ 



where the matrices 


ofM 




/II 


1 


1 


1 


0\ 




1 


1 


1 


1 







1 


1 


1 


1 


1 


Ai = 


1 


1 


1 


1 


1 







1 


1 


1 


1 







1 


1 


1 


1 1 




\0 


1 


1 


1 


11/ 




/II 


1 


1 





0\ 




1 1 


1 


1 


1 







1 1 


1 


1 


1 





A3 = 


1 


1 


1 


1 







1 


1 


1 


1 


1 




1 


1 


1 


1 


1 




\0 


1 


1 


1 


10/ 



In [PJBj the values p{A4) and p{A4) were computed approximately to the following 
accuracy: 

4.525 < p{M) < 4.6105 ; 4.72 < p{M) 

Algorithms (P) and (L) determine their precise values: 

,1/2 



< 4.8. 



p{M) = [p{AiA2)y 
MM) = [p{A2A3)] 



1/2 



4.61047781. 
4.72204513. 



The JSR computation. Algorithm (L) starting with the product 77 = A2 A3 
terminates after 4 steps producing the extremal infinite polytope 7-3 = co-{vi}}li, 
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where vi is the leading eigenvector of 77, V2 = A3V1 (step zero); V3 = Aivi, 
V4, = A2V1, V5 = A1V2, and ve = A3V2 (first step); V7 = A^Vi, vs = Aive, 
vg = A2Ve, and vio = A3V6 (second step) ; vn = Aivs and V12 = A2VS (third 
step). See the corresponding cyclic tree in Figure [9l Thus 11 = A2 A3 is an s.m.p. 




Fig. 10 The vertices {vi}j2^-^ for LSR computation of the ternary Eulor partition function 

The LSR computation. Algorithm (P) starting with the product 77 = Ai A2 
terminates after 4 steps, the extremal polytope P3 = co-^.{vi}l2=l■ The vertex vi 
is the leading eigenvector of 77, V2 = A2V1 (step zero); V3 = Aivi, V4, = ^3^1, 

= A2V2, and VQ = A3V2 (first step); V7 = A1V3, vs = A2V3, vg = ^31^3, 
viQ = A1V4, vii = A2V4, V12 = A1V5, and 1113 = ^21^5 (second step); V14 = 
Aivg, vi5 = ^21^9, and viq = A3Vg (third step). See the corresponding cyclic 
tree in Figure [TOl Therefore 77 = Ai A2 is an s.l.p. 

9 Numerical results for randomly generated matrices 

In this section we report some results obtained for families consisting of a pair of 
random matrices of variable dimensions d. 

The results show that the computation complexity increases significantly as the 
dimension increases but also confirm the effectiveness of the method for computing 
the joint and the lower spectral radius of nonnegative matrices. We expect in 
general that the reachable dimension for a computation in a reasonable time might 
be quite high for a set of operators sharing an invariant cone. 

First we consider the general case of two random matrices with normally dis- 
tributed entries. The generated random matrices are scaled to have equal spectral 
norm. This aims to reduce the number of cases were the s.m.p is the matrix with 
larger spectral radius. The first column of Table [2] gives the dimension, the second 
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JSR 






JSR 




J 

a 


-U- ? + 

# Its 


^ vertices 


s.m.p. 


a 


# its 


^ vertices 


S.m.p. 


5 


3 


14 


A A 


6 


4 


26 




5 


7 


23 




6 


9 


51 


A A 

A1A2 


r 





QV 
O ( 


A 

A\ 





5 


38 


a2 a 


7 


17 


100 


A\ 


8 


19 


117 


a3 A aA A 
A-{A2A\A2 


7 


12 


140 


a'^ A A A 

A\A-2A\A-2, 


8 


8 


49 


Ai 


/ 


0/1 

Z4 




A'i A'i 


Q 


12 


75 


A zl 3 


9 


18 


177 


zl8 A 
A\A2 


10 


16 


239 


AlA^ 


n 

y 


1 Q 

lo 


1 vo 
1 i z 


zl3 AAA 




9 


109 


A 

Ai 


9 


10 


129 


A-, 


10 


24 


408 


{A\A2fA2 


11 


20 


707 


A\Al 


12 


31 


1539 


A1A2AIAI 


11 


14 


340 


AIA2A1A2 


12 


9 


211 


A1A2 


11 


12 


183 


A\A2 


12 


13 


215 


MAl 


15 


18 


715 


A\A2AxA\ 


20 


21 


1539 


A1A2 


15 


14 


570 


A\A2 


20 


16 


1219 


AiAl 


15 


14 


390 


A2 


20 


16 


1247 


A^A2 



Table 2 Computation of the JSR for random pairs of matrices with equal norm. 



column the number of iterations for Algorirhm (R) to converge, the third column 
provides the number of vertices of the extremal polytope and the last column gives 
the correspondent s.m.p.; we immediately observe that the complexity (in terms 
of iterations and number of vertices) rapidly increases with the dimension. Dealing 
with two 20 X 20 matrices can be considered a challenging computational problem. 

Then we consider in Table [3] randomly generated nonnegative matrices still 
scaled to have the same norm. 







JSR 






LSR 




d 


# its 


# vertices 


s.m.p. 


# its 


# vertices 


s.l.p. 


10 


3 


6 


A1A2 


4 


6 


AiAl 


10 


3 


4 




4 


5 


A2 


10 


4 


6 


A1A2 


7 


15 


A\A% 


10 


6 


11 


AlAl 


3 


6 


A1A2 


10 


4 


8 


AiAl 


5 


9 


A\A2 


20 


4 


7 


A2 


4 


6 


Ai 


20 


4 


6 


A1A2 


5 


9 


AiAl 


20 


6 


14 


AIA2 


3 


4 


A1A2 


20 


5 


11 


AiAl 


6 


14 


A\A2 


20 


5 


9 


A1A2 


3 


4 





Table 3 Computation of the JSR and of the LSR for random nonnegative pairs of matrices. 



Finally we consider binary matrices and vary the density of the number of zero 
entries. We scale the pairs of matrices to have the same spectral radius; note that 
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in some cases either the s.m.p. or the s.l.p. are the starting matrices so that there is 
no guarantee of the convergence of the algorithm we propose. Nevertheless, when 
we report the number of iterations and of vertices we imply that the algorithm 
has converged in a finite number of steps. 

Tables Eland [5] report the results obtained for pairs of matrices respectively of 
dimension d = 50 and d = 100. Whenever both Ai and A2 are either s.m.p. 's or 
s.l.p. 's we indicate both the numbers of iterations/ vertices taking either Ai or A2 
as optimal product. 







JSR 






LSR 




density 


# its 


# vertices 


s.m.p. 


# its 


# vertices 


s.l.p. 


0.2 


9 


55 


AlAl 


4 


8 


Ai and A2 


0.2 


5 


17 


A1A2 


5 


10 


AIA2 


0.2 


8 


24 


AiAl 


4(4) 


6(6) 


Ai and A2 


0.2 


5 


16 


AIA2 


4(5) 


6(8) 


Ai and A2 


0.2 


14 


59 


AlAl 


5 


10 


A1A2 


0.5 


4 


8 


A1A2 


4 


10 


A1A2 


0.5 


5 


14 


AIA2 


4(3) 


5 (4) 


Ai and A2 


0.5 


6 


15 


AlAl 


6 


17 


AlAl 


0.5 


5 


16 


A1A2 


4(4) 


6(5) 


Ai and A2 


0.5 


6 


20 


AIA2 


5 


9 


A1A2 


0.75 


5 


16 


AlAl 


5 (7) 


12 (14) 


Ai and A2 


0.75 


4 


8 


A1A2 


6 


16 


A1A2 


0.75 


5 


11 


AlAl 


6 


19 


AlAl 


0.75 


5 


16 


AlAl 


11 


170 


AiA% 


0.75 


5 


12 


A\A2 


5 (6) 


13 (12) 


Ai and A2 


0.9 


4(5) 


8(9) 


Ai and A2 


4 


8 


A1A2 


0.9 


5 


9 


A\A2 


6 


4 


A1A2 


0.9 


3 


4 


A1A2 


7(8) 


11 (12) 


Ai and A2 


0.9 


4 


11 


AlAl 


4 


7 


A1A2 


0.9 


7 


14 


AlAl 


8 (8) 


13 (11) 


Ai and A2 



Table 4 Computation of the JSR and of the LSR for random pairs of binary matrices of 
dimension d = 50. 



Some comments are necessary. The computations for the general case with 
d = 10 need usually a few minutes. The computations for the general case with 
d = 20 need usually between half an hour and one hour of computation but for 
some examples till 8 hours (in a standard laptop with 15 processor) . We hypothesize 
that the overall computation depend on several factors, not only the length of the 
spectrum maximizing product but also the ratio between the leading eigenvalues of 
the family and the closer ones that is eigenvalues of products which have modulus 
close to 1 and on the distribution of the vertices of the extremal polytope in all 
the orthants. 

In the nonnegative case all the vertices lie in the nonnegative orthant and this 
determine a much lower complexity. The presence of quasi-optimal products that 
is products with eigenvalues close to 1 is a factor of slowdown also in this case. 
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# Its 
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s.m.p. 
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7^ vertices 


s.l.p. 


0.2 


6 


24 


AfA2 


6 


31 


AiA^ 


0.2 


6 


23 


A1A2 


6 


28 


a2 a2 


n o 


i 


ov 
z ( 


A A3 





on 


A A 
A1A2 


n o 





01 

zi 


A 42 

AlAy 


/ 


z4 


A2 a 
A-f A2 


0.5 


5 


10 


A A 

A\A2 


5 


15 


A a2 
A1A2 


0.5 


6 


17 


AIA2 


4 


8 


A1A2 


0.5 


6 


18 


AlAl 


5 


16 


AIA2 


0.5 


6 


22 


MAl 


4(6) 


9 (14) 


Ai and A2 


0.8 


4 


7 


AxA2 


4 


7 




0.8 


7 


18 


AIA2 


6 


14 


AlAl 


0.8 


5 


14 


AiAl 


9(7) 


14 (16) 


A\ and A2 


0.8 


5 


12 


AIA2 


5 


12 


AlAl 



Table 5 Computation of the JSR and of the LSR for random pairs of binary matrices of 
dimension d = 100. 



For binary matrices we observe from the experiments that the behavior of the 
algorithm slightly depends on the density of the zero entries and also on the dimen- 
sion. This implies that we are able to compute the joint spectral characteristics of 
possibly large binary matrices with any density. 

In some applications the families of the matrices have the same spectral radii 
and it may happen that the optimal products are exactly the matrices themselves. 
In our experience there are cases where we have been able to compute an extremal 
polytope invariant set by starting from any optimal matrix of the family (see also 
the tables of results) but we have also encountered cases where the algorithm 
has not terminated finitely. We think that an interesting open problem is that 
of balancing leading eigenvectors associated to different products (which are not 
powers or cyclic permutation one of the other). 

Although it is true that this situation is not generic there are some applications 
where it naturally occurs. We leave this topic to a future investigation. 

10 Appendix 

We give here the proofs of two main results of this paper (Theorems [3] and [7]) and 
details about the proof of Theorem [8l 



10.1 Proof of Theorem m 

We give the proof for Algorithm (R), the proofs for Algorithms (C) and (P) are 
analogous. We use two auxiliary results. 

Lemma 12 Let us have a cyclic tree T with a root B generated by an irreducible 
word b = dn-.-di; then for any word a, which is not a power of b, we have 
a" ^ B. 
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Proof . Let I be the length of a and p be and the greatest common divisor of I 
and n. If Z = kn for some integer k, then the words a and b'° must have different 
letters at some position, otherwise a = b*^. Therefore a ^ B, and so a" ^ B. 

If I is not divisible by n, then p < n, and there exists an index j such that 
dj+p / dj, otherwise b is a power of the word dp . . . di, which contradicts the 
irreducibility. The Diophantine equation Ix — ny = p has a solution {x,y) £ 
such that < x < 71 — 1. Since the words a^+^ and b^+^ have different letters at 
the position Ix + j, we have a"^"*"^ ^ B, and hence a" ^ B, because n> x + 1. 

□ 

Lemma 13 LetT he the cyclic tree generated by the product U. If Algorithm (R) 
terminates within k steps, then there is e > such that all nodes of T of level > k 
are in the polytope (1 — e) Pk-i- 

Proof. If the algorithm terminates after k steps, then AjPk-i C Pk-i for all 
j = 1, . . . ,m. Moreover we have Uk = 0, which means that every node v (z T of the 
level k belongs to a dead branch, and therefore v G int Pk-i- The total number of 
nodes of level k is finite, hence all of them are in (1— e) Pk-i for some e > 0. If w is a 
node of a bigger level r > k, then v = Rvq, where R G A4 and vq is some node 
of the fcth level. Since uq € (1 — e) -Pjt-i, we have v G (1 — e) RPk-i C {1 — s) Pk-i, 
because RPk-i C Pk-i- 

□ 

Proof of Theorem 15] 

Necessity. Consider the cyclic tree T generated by the product 77. Assume the 
algorithm terminates after k steps. By Lemma [13] all nodes of levels at least k 
belong to (1 — e) Pk-i, where e > is fixed. For every product C, which is not a 
power of 77i, the node C"ui does not belong to the root (Lemma I12p. Hence for 
each Vi from the root and for every product C G M' that is not a power of a cyclic 
permutation of 77, the level of the node C^^'^Vi is bigger than k. If v is not in 
the root, then this level is bigger than In + k > k. Thus, C"^'^w G (1 — e) T-'^-i 
for each node f G T, and hence, for each vertex v of the polytope Pk-i- This 
yields that C+'^Pk-i C (1 - e)Pk-i. Therefore, p(C"+'=) < 1 - e, and so 
p(C) < (l — e)^^^^^''\ Consequently 77 is dominant. 

Let us now show that 1 is its unique and simple leading eigenvalue. Since for 
i 7^ 1 the product 77i is not a power of 77i, it follows that the node UiVi does not 
belong to the root fLemmall3|l. Hence, the level of the node Il^'^^Vi is bigger than 
k. If V is not in the root, then the level of U^'^'^v is bigger than k as well. Thus, 
jj^+f'y g (2. — e) Pk-1 for all vertices v of T^-i, except for v = ±wi. For any 
eigenvector it ^ ui of the operator 77i take the one-dimensional subspace [/ C M*^ 
spanned by u (if u is complex, then U is the two-dimensional subspace spanned by u 
and by its conjugate). Since vi ^ U, it follows that ni{Pk-iC\U) C int {Pk-iC\U), 
where the interior is taken in U . This implies that the spectral radius of 77i on 
the subspace U is smaller than 1. Thus, all eigenvalues of 77i different from 1 are 
smaller than 1 by modulo, and the eigenvalue 1 has a unique eigenvector. Hence, 
the leading eigenvalue 1 is unique and has only one Jordan block. The dimension of 
this block cannot exceed one, otherwise ||77{°|| — > oo as fc — > oo, which contradicts 
the nondefectivity of the family Ad. Therefore, the eigenvalue 1 is simple. 
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Sufficiency. The proof uses similar arguments as this given for the Small CPE 
Theorem in |GWZj . to which we refer the reader. 

Assume U is dominant and its leading eigenvalue is unique and simple. If the 
algorithm does not terminate, then the tree T has an infinite path of ahve leaves 
v^'^^ — > D^^' — > . . . (the node d^*' is on the ith level) starting at a node = Vp 
from the root. For every r we have w*''^ ^ int Pr-i- Hence u'-''-' ^ int Pk for aU k < r. 
Since the family M is irreducible, it follows that G int P], for some fc, and hence 
the polytope Pk defines a norm || • \\k in . For this norm > 1 for all r > k. 

On the other hand, A4 is nondefective, hence the sequence {v^^^} is bounded. Thus, 
there is a subsequence ri > fc, that converges to some point v £ R'^. 

Clearly, ||u||fc > 1. For every i we have i;''''^ = RiV^^^'' andu^'"'+^^ = Ci where 
Ri e M'''~''^ , Ci e M'''+^~''\ Denote by ci[M] the closure of the semigroup of 
all products of operators from A4. Since this semigroup is bounded, after possible 
passage to ^ subsequence, it may be assumed that Ri and d converge to some 
R,C £ c£[M] respectively as i — > oo. We have Cv = v, hence p{C) > 1, which, 
by the domination assumption, implies that there is j G {1, . . . ,n} such that C 
belongs to c^[/7j], which is the closure of the semigroup {(77j)'}ggN- Moreover, 
since the leading eigenvalue of 77 is unique and simple, we see that v = \vj, where 
A G m. We have, \\vj\\k = 1 and \\v\\k > 1, hence |A| > 1. Thus, Rv'-'"''^ = Xvj. The 
nodes vj and v^^'' = Vp are both from the root, hence there is a product S such that 
Svj = v^'^\ Taking into account that v^''^^ = Riv^'^\ we obtain RiSRv^^'^ = Xv^'^\ 
Hence p{RiSR) > |A|, and we conclude that A = ±1 and that RiSR G c^[77i] 
for some i. This yields t;^^-' = fivi, where Vi is the corresponding vector from the 
root, /.t G K. Since \\vi\\ = 1 and Hf''"'^'' || > 1, we have \fi\ > 1. The elements Vi 
and v^^^ = Vp are both from the root, hence Ap ■ ■ ■ AiVi = Vp, and consequently 
AsAp • • • Ai Vi = ' for some As e M. Note that dg ^ dp+i, because the node 
v^^'^ is not in the root. Therefore, the product Q = AsAp ■ ■ ■ Ai does not coincide 
with 77i, and its length is at most n. Thus, Qvi= fiVi, hence p{Q) = 1 and /x = ±1. 
Thus, Q has spectral radius 1 and the leading eigenvector Vi, therefore Q G c£[77i]. 
On the other hand the length of Q does not exceed n, hence Q = Hi, which is a 
contradiction. Hence, the algorithm terminates within finitely many steps. 

□ 



10.2 Proof of Theorem □ 

We use several auxiliary results. The proof of the following lemma is similar to the 
proof of Lemma 1131 

Lemma 14 Let T be the cyclic tree generated by the product 11 . If the algorithm 
terminates within k steps, then there is s > such that all vertices of the tree of 
level > k belong to the infinite polytope (1 + e) Pk-i- 

Lemma 15 jV^ If an operator B has an invariant cone K , then for every its 
eigenvector from int Ti' the corresponding eigenvalue equals to p{B). 

Proof of Theorem IT] 
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Necessity. Consider the cyclic tree T generated by the product 77. If the al- 
gorithm terminates after k steps, then by Lemma [T4l all vertices of levels at least 
k belong to (1 + e) Pk-i- For an arbitrary product C, which is not a power 
of 77i, the node C^Vi does not belong to the root fLemma I12|) . Hence for ev- 
ery Vi from the root the level of the node C'^'^^Vi is bigger than k. If v is 
not in the root, then the level of is bigger than k as well, consequently 

(j'^+^y G (l + £)7-fe_i for each node v & T, and hence for every vertex of Pfc. 
This yields that C+'Tk-i C {1 + e) Pk-i, therefore p(C"+'=) > 1 + e, and so 
p(C) > (l + e)^^''"^'^^ This holds for every product C that is not a power of 77 
or of its cyclic permutations, which completes the proof. 

Sufficiency. Assume the converse: the product 77 is under-dominant, but the 
algorithm does not produce an extremal infinite polytope. This means that the 
tree T has an infinite path of alive leaves 

„(o) ^ „{i) ^ 

. starting at a vertex 

from the root. Since the family M is eventually positive, it follows from 
Lemma |9] that there exists an internal invariant cone K, which, moreover, con- 
tains all leading eigenvectors of products of operators from A4. Hence, K contains 
the root of T, and therefore, it contains all the nodes v^''\ For every r we have 
t;^''^ ^ intPr-i. Hence t)*-'"-' ^ intPfc for all k < r. Let </(•) be the antinorm defined 
by the infinite polytope Pk: g{x) = sup •[ A | X~^x G Pk }. Since a concave func- 
tion is continuous at every interior point of its domain (see, for instance, |MT] ) . 
it follows that g is equivalent to every norm and to every antinorm on the interior 
cone K. In particular, there are positive constants ci, C2 such that 

ci/(x) < g{x) < C2f{x) , X e K 

where / is an invariant antinorm for A4 (Theorem [6]) . For arbitrary r we have 
p(t;('')) < 1. On the other hand, since / is invariant and f{vk) = 1 for all k, 
we have f{v) > 1 for every node v of the tree. In particular, f{v^^^) > 1. 
Thus, ci < g(w(''^) < 1 for aU r. Since g is equivalent to each norm on K, 
we see that the sequence {v^^''} is bounded, and hence there is a subsequence 
, ri > 1 , that converges to some point v G K"^. Clearly, v £ K and 
ci < giv) < 1. For every i we have w'-'''-' = 7iiW*-'''' and = CiV^^^\ where 

Ri e M , Ci e M The sequence {v'^'''^} is contained in K, bounded, 

and separated from zero, hence by Lemma 1111 the sequences of operators {Ri} 
and {Ci} are both bounded. Therefore, after a passage to subsequences it may be 
assumed that these two sequences converge to some R,C £ c£[AI] respectively as 
z — > cxD (see the proof of Theorem [3] for the definition of ci?[A^] and c£[77]). We 
have Cv = v. Since v G intRiJ., if follows from Lemma [131 that v is the leading 
eigenvector of C. Consequently, p{C) = 1, which, by the domination assumption, 
implies C G c^[77j] and v = Xvj for some j = 1, . . . ,n, and A > 0. Since g{vj) = 1 
and g{v) < 1, it follows that A < 1. Thus, Rv'''^'-^ = Xvj. The elements vj and 
v^^^ = Vp are both from the root, hence there is a product S such that Svj = v^^\ 
Taking into account that v'^^^^ = Riv''^\ we obtain RiSRv^^^ = Xv^^'. Again 
invoking Lemma [TS] we conclude that is the leading eigenvector of RiSR. 
Hence p{RiSR) = A, and we see that A = 1, hence RiSR G c^[77i] for some i. 
This yields v^^^ = fivi, where Vi is the corresponding vector from the root, fi > 0. 
Since g{vi) = 1 and g{v^^^) < 1, we have n < I. Elements Vi and = Vp are 
both from the root, hence Ap ■ ■ ■ AiVi = Vp, and consequently AgAp ■ ■ ■ AiVi = w'^-' 
for some As G A4. Note that ds ^ dp^i, because the vertex v''^^ is not in the 
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root. Therefore, the product Q = AsAp ■ ■ ■ Ai does not coincide with 77i, and its 
length is at most n. Thus, Qvi = fivi, and by Lemma llSI o(0) = fi. Consequently, 
/i = 1 and Vi is the leading eigenvector of the operator Q G c£[i7i]. On the other 
hand the length of the product Q does not exceed n, therefore Q = Hi, which is 
a contradiction. 

□ 



10.3 The 20 x 20-matrices Ai,A2 for the problem of overlap-free words of i)8.1l 
and the proof of Theorem [51 

We write the two matrices Ai, A2, associated to the problem discussed in H8A 
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respectively. 

To give the rigorous proof of the theorem it now suffices to list all the vertices 
of the extremal polytopes obtained by applying Algorithms (P) and (L). 

Proof of Theorem [H Denote 11 = A1A2. To show that 



p(Ai,A2) = [p(77) 



1/2 



it suffices to present an extremal polytope P for the operators Ai = [p(-/7)] Ai 

and A2 = [p(77)]"^^^A2. This polytope is _P = CO- where the first 

vertex vi is the leading eigenvector of 77, and the other vertices are 



1)2 = A2V-i, V3 = 

VS = Aivs , vg = 

fl4 = AlVQ, Vl5 

f20 = A2Vg, V21 

f26 = AlVie, V27 

f32 = A2V14, D33 

V3S = A2V22, V3g 

1)44 = ylill33, Vis 
VBO = ^1^43, 



AiVi, V4 = 

Aivs , f 10 

= Aivio, Via 

- A2V10, V22 

- A1V17, V2S 

- A2V1Q, 1134 
= -4.1 ''as, ''40 

= ^^-"25, V46 

= A1V45, V52 



A1V2, fs = 

= A2V4,, V\l ■ 

- Mvw, 1)17 ■ 
= A2V\\, V23 ■ 

- Aivig, V29 ■ 

= A2V17, D35 : 

= A\V2&, D41 ■ 

- A1V46, D53 : 



A2V2, V(j = 

= A2V5, 1)12 ■■ 

= A\vi2, D18 ■ 

= A2V12 , 1)24 ■ 

= A1V21, 1)30 ■ 

= A2Vig, 1)36 ■ 

= A\V27, D42 ■ 

= ^2'«30> 1^48 ■ 

= A2D39, D54 ■ 



A2VZ, V7 = A1V4, 

= A2D6, D13 = A1D8, 

= A2V7, D19 = A2D8, 

= A1V13, D25 = A1D14, 

= A1D22, D31 = A2D13, 

= A2V20, D37 = A2V21, 

= A1D30, D43 = A1D32, 

= A2D33, D49 = A1D39, 
= ^21)45- 
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Now let n = AiAl". To prove that 

p{Ai,A^) = [pin)]'/'' 

it suffices to present an extremal infinite polytope P for the operators Ai = 
[p(n)]~'/"Ai and I2 = [pin)]~'/"A2. This polytope is P = co{{vi}jZ\) + 
Kfi, where the first vertex vi is the leading eigenvector of 11 , and the other vertices 
are 



V2 — 


- A2V\ , 


V3 - 


— A2V2 5 


V4 = 


= A2V3, 


"5 - 


= ^2"4, 


"6 = 


= ^2"5, 


V7 


- yl2"6. 


V8 = 


= ^42117, 


V9 - 


= A2"8, 


"10 


= A2V9, 


"11 


= A2"10, 


"12 


= AlVl, 


ri.i 


= A1V2, 


V14, 


= ^1^3, 


fl5 


= A\V4, 


"16 


= Awe, 


"17 


= A1V7, 


"18 


= Aivs, 


"19 


= A1V9, 


V20 


= Mvu), 


V21 


= M'uii, 


"22 


= A1V12, 


"23 


= Alvri, 


"24 


= Aivis, 


"25 


= 41 "19, 


V26 


= MV20, 


V27 


= Aru2i, 


"28 


= A2V12, 


"29 


= 42"18, 


'"30 


= Mvig, 


"31 


= 42 "20, 


V32 


= A2V21, 


"33 


= A1V22, 


"34 


= 41 "27, 


"35 


= 4,1 ''so. 


"36 


= AIV31, 


"37 


= 41 "32, 


V38 


= MV27, 


"39 


= 42 "30, 


"40 


= 42"31, 


"41 


= A2V32, 


"42 


= AIV34, 


"43 


= 41 "37, 


V44 


= MV38, 


"45 


= A1V40, 


"46 


= A1V41 , 


"47 


= 42 "34, 


"48 


= A2V3S, 


"49 


= 42 "40, 


f50 


= A1V42, 


"51 


= A1V43, 


"52 


= A1V44, 


"53 


= AlV48, 


"54 


= A2V42, 


"55 


= A2V43, 


V5(i 


= A2V44, 


■"57 


= A2V47, 


"58 


= A2V4S, 


"59 


= ^1"50, 


"60 


= AiVo2, 


"61 


= A1V54, 


f62 


= A2V50, 


«63 


= A2V52, 


"64 


= A2V57, 


"65 


= A1V59, 


"66 


= Alvao, 


"67 


= ^1"62, 


V68 


= A2V59, 


"69 


= A2"60, 


"70 


= A2VQ2, 


"71 


= AlV65, 


"72 


= 41 "67, 


"73 


= A1V68, 


V74 


= AlV69, 


"75 


= Aii;70, 


"76 


= A2V65, 


"77 


= A2V67, 


"78 


= A2"68, 


"79 


= Aij;7i, 


"80 


= MV72, 


"81 


= 41 ■"73, 


"82 


= AiV76, 


"83 


= A1V77, 


"84 


= 41 ''78, 


"85 


= 42 "71, 


V86 


= A2V72, 


"87 


= A2V73, 


"88 


= A2V7e, 


"89 


= A2V7S, 


"90 


= 41 "79, 


"91 


= 41 "81, 


V92 


= AiV82, 


"93 


= 41 "84, 


"94 


= ^1"85, 


"95 


= AiVs7, 


"96 


= ^1^88, 


"97 


= -42"79, 


V98 


= A2V8I , 


"99 


= ^2 "82, 


"100 = ^2 "85, 


"101 


. = -42"87, 


"102 = ^2"88, 


"103 


= ^1"92 


V104, 


= ^2^100 ■ 























The proof is completed by routine computations. 

□ 
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